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The  accuracy  of  modified  Sachs  and  Ledsham-Pike  scaling  of 
peak  shock  hydrodynamic  variables  from  a  nuclear  burst  in  air  is 
evaluated.  The  modified  Sachs  and  Ledsham-Pike  methods  are 
corrections  applied  to  the  similarity  transform  which  is  used  to 
compute  shock  overpressures  and  related  variables  for  infinite 
homogeneous  atmospheric  ambient  conditions.  This  similarity 
transform  no  longer  applies  when  the  burst  and  target  are  located 
at  different  altitudes,  and  the  modified  Sachs  and  Ledsham-Pike 
corrections  are  applied  to  account  for  these  varying  ambient  con¬ 
ditions.  In  this  dissertation,  the  modified  Sachs  shock  posi¬ 
tions,  dynamic  pressures,  and  overpressures  and  the  Ledsham-Pike 
overpressures  are  compared  to  a  fully  two  dimensional  flux- 
corrected  transport  finite  difference  solution.  The  flux- 
corrected  transport  method  maintains  a  sharp  and  steady  shock 
with  no  oscillations.  The  numerical  calculations  include  a  real 
atmospheric  model,  the  Doan  and  Nickel  equation  of  state  of  air, 
and  radiation  energy  losses.  It  is  found  that  the  modified  Sachs 
and  Ledsham-Pike  overpressures  are  essentially  identical,  and 
that  their  overpressures  compare  with  the  numerical  results  more 
favorably  in  the  ascending  directions  than  in  the  descending 
directions.  A  new  L'  isham-Pike  correction  factor  is  presented 
for  the  descend  directions.^ 
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Chapter  I.  Introduction 


(•>* 

When  a  nuclear  weapon  is  detonated  in  the  atmosphere  the 
immediate  effect  is  the  liberation  of  a  large  amount  of  energy  in 
a  small  volume.  Between  seventy  and  eighty  percent  of  the 
weapon's  yield  appears  in  the  form  of  X-rays  (see  Glasstone  and 
Dolan,  1977) ,  and  since  the  mean  free  path  of  these  X-rays  in  air 
is  a  few  meters,  they  are  absorbed  within  several  meters  of  the 
burst  point.  The  air  becomes  ionized,  and  the  resulting  plasma, 
which  includes  bomb  debris,  is  called  the  fireball.  Initially, 
the  radiation  energy  from  the  weapon  is  greater  than  that  re¬ 
quired  to  completely  ionize  the  fireball's  constituent  atoms. 
The  X-rays  then  stream  through  the  plasma  and  are  absorbed  at  the 
fireball's  boundary,  ionizing  successive  layers  of  air  near  the 
boundary.  The  fireball  thus  grows  until  the  radiation  energy 
from  the  weapon  can  no  longer  reach  the  surface  of  the  fireball. 
This  point  is  called  "burn-out".  The  fireball  continues  to  grow 
diffusively  after  burn-out  by  X-ray  re-emission  and  re¬ 
absorption,  but  the  growth  rate  continuously  slows  down.  At  some 
time  during  the  diffusive  growth  phase  of  the  fireball,  the  po¬ 
tential  shock  wave  speed  from  the  blast  is  greater  than  the  rate 
of  growth  of  the  fireball.  At  this  time  the  shock  wave  emerges 
from  the  fireball,  and  this  phenomenon  is  called  "hydrodynamic 
separation" . 

After  hydrodynamic  separation,  the  shock  wave  propagates 
outward,  and  its  motion,  as  well  as  the  motion  of  the  fluid 
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behind  the  '•hock  front,  is  governed  by  the  partial  differential 
equations  of  hydrodynamics.  These  are  continuity  equations  ex¬ 
pressing  the  conservation  of  mass,  momentum,  and  energy  (see 
Zel'dovich  and  Raizer,  1966;  or  Harlow  and  Amsden,  1971),  and  in 
one  dimension  they  can  be  written 

&  *  h  £rV  +  S  =  o 


where 


f 

o=  n 

E 


(1-2 


pr 

f(0)=  fVl 

E-V 


s  = 


(1-4 


and  v  =0,  1,  2  for  Cartesian,  cylindrical  and  spherical  geometry 
respectively.  For  a  complete  evaluation  of  the  hydrodynamics 
variables  as  functions  of  time  and  space,  the  three  equations 
represented  by  (1-1),  together  with  an  appropriate  equation  of 


state,  must  be  simultaneously  solved  for  the  mass  density 


momentum  density 


v,  energy  density  E,  and  nr es sure  p.  In  the 


most  general  case,  the  complete  multidimensional  hydrodynamics 
equations  would  be  solved  with  initial  conditions  appropriate  for 
each  atmospheric  condition  of  interest.  These  extremely  expen¬ 
sive  repeated  calculations  can  be  (and  are)  avoided  by  a  scaling 
method  which  uses  a  single  complete  solution  of  the  one¬ 
dimensional  equations  ( I — 1 )  for  a  given  yield  in  a  uniform  homo¬ 
geneous  atmosphere,  and  then  transforming  this  solution  to  anoth¬ 
er  yield  and  atmosphere.  The  similarity  transformation  which 
allows  this  procedure  is  known  as  Sachs  scaling,  and  it  is  used 
extensively  by  the  nuclear  survivability  and  vulnerability 
analysts.  The  "given"  complete  solution  is  usually  for  a  one- 
kiloton  burst  in  a  standard  temperature  and  pressure,  sea  level 
(but  no  ground  effects)  atmosphere.  Since  this  solution  is  wide¬ 
ly  used  in  the  nuclear  effects  community,  the  Air  Force  Weapons 
Laboratory  (AFWL)  has  prepared  a  compendium  of  one-dimensional 
numerical  solutions  and  experimental  measurements,  all  scaled  to 
the  one  kiloton  yield  and  standard  sea  level  conditions  (Needham 
and  Havens,  1975).  This  averaged  survey  of  results  has  become 
known  as  the  AFWL  1  kt  standard. 


With  the  AFWL  1  kt  standard  as  the  known  solution,  the  Sachs 
scaling  laws  can  be  used  to  transform  all  the  flow  variables  to 
any  other  yield  and  homogeneous  atmosphere;  however,  the  similar¬ 
ity  transformation  gives  no  information  about  the  case  when  the 
burst  is  at  one  altitude  and  the  target  is  at  another  altitude. 
In  this  case,  the  shock  travels  through  a  variable  density  and 


pressure  atmosphere  and  these  spatially  variable  external  condi¬ 
tions  have  a  pronounced  effect  on  the  shock  values  over  inter¬ 
mediate  and  large  distances  (several  kilometers).  To  account  for 
these  effects,  a  procedure  known  as  modified  Sachs  scaling  is 
used.  This  procedure  is  defined  in  the  next  section  and  has  no 
theoretical  basis,  but  when  compared  to  measured  data,  it  has 
been  found  to  give  "reasonably  reliable"  answers  (Moulton,  1960). 
A  second  correction  procedure  is  available,  but  has  received 
little  attention.  It  is  the  Ledsham-Pike  correction,  and  is  used 
to  correct  overpressures  values  only  (Ledsham  and  Pike,  1951). 
The  Ledsham-Pike  correction  factor  has  a  theoretical  basis,  but 
its  values  must  be  experimentally  measured  or  computed  from  a 
fully  two-dimensional  solution  of  the  equations  of  hydrodynamics. 

This  dissertation  is  a  critical  evaluation  of  the  modified 
Sachs  and  Ledsham-Pike  correction  techniques.  They  are  evaluated 
for  one  set  of  initial  and  external  conditions  by  comparing  their 
results  to  a  two-dimensional  finite  difference  solution  which 
uses  the  method  of  flux-corrected  transport.  The  particular  case 
chosen  for  comparison  was  a  one  megaton  burst  at  15  kilometers 
altitude.  The  next  section  is  a  discussion  of  the  Sachs  and 
modified  Sachs  scaling  methods  and  the  Ledsham-Pike  correction 
factor.  It  is  followed  by  a  description  of  the  numerical  calcu¬ 
lations  required  for  the  comparison  including  a  discussion  of 
finite  differences  and  flux-corrected  transport.  The  fourth  sec¬ 
tion  compares  the  scaling  law  predictions  to  the  numerical 
results,  and  this  is  followed  by  the  conclusions. 


Chapter  II.  Scaling  Laws 


To  predict  shock  densities,  velocities  and  pressures  for 
different  initial  conditions  (yields)  and  ambient  or  atmospheric 
conditions  requires  a  solution  of  the  partial  differential  equa¬ 
tions  of  hydrodynamics  or  experimental  measurement  of  the  proper¬ 
ties  for  each  set  of  conditions.  Fortunately,  we  can  use  simi¬ 
larity  theory  to  develop  scaling  relationships  to  avoid  such  a 
costly  procedure.  The  fundamental  principle  underlying  these 
relationships  is  that  the  same  conservation  laws  govern  the  fluid 
flow  and  shock  wave  mechanics  for  all  explosions.  For  distances 
far  enough  removed  from  the  explosion  that  the  mechanism  by  which 
the  explosion  energy  is  transfered  to  the  air  is  not  important, 
the  only  governing  parameters  are  total  explosion  energy  and  the 
ambient  conditions  of  the  air.  Thus,  explosions  in  air  obey 
principles  of  similitude: 

"Two  phenomena  are  similar,  if  the  characteristics  of 
one  can  be  obtained  from  the  assigned  characterictics 
of  the  other  by  a  simple  conversion,  which  is  analogous 
to  the  transformation  from  one  system  of  units  of 
measurement  to  another."  (Sedov,  1959) 

The  transformation  functions,  or  similarity  ratios,  are  called 
scaling  factors,  and  are  generally  derived  from  the  hydrodynamics 
equations  when  they  are  reformulated  as  relations  between 


nondimens ional  quantities.  (Dimensional  quantities  depend  on  the 
units  of  measurement  used  in  their  numerical  definition  -  nondi- 
mensional  quantities  do  not  depend  on  the  measurement  scale  and 
their  numerical  values  are  the  same  in  all  systems.) 

One  of  the  earliest  scaling  laws  for  high  explosive  detona¬ 
tions  was  proposed  by  Hilliar  in  1919  (see  Moulton  et  al,  1971)  . 
Since  the  energy  available  in  a  charge  is  directly  proportional 
to  the  mass,  and  the  cube  root  of  the  mass  is  directly  propor¬ 
tional  to  the  characteristic  length  of  the  volume  containing  the 
mass,  Hilliar  concluded  that  the  scaling  factors  for  radius  and 
time  for  two  similar  explosions  of  different  yields  but  in  iden¬ 
tical  air  media  must  be  the  cube  root  of  the  ratio  of  the  two 
explosives'  masses  (or  energies).  Thus,  for  example,  the  peak 
pressure  p2  of  a  TNT  explosion  of  weight  w2  at  radius  r2  is 

Pjfti)  -  k£r>)  <ii-i 

where 


and  p  is  the  peak  pressure  of  a  "reference"  explosion. 


Sachs  and  Modified  Sachs  Scaling 

Hilliar's  yield  scaling  laws  developed  for  high  explosive 
shock  waves  apply  to  nuclear  explosions  in  an  infinite  homo¬ 
geneous  atmosphere  when  the  weight  of  the  nuclear  explosion  is 
defined  by  the  "TNT  equivalent"  (or  when  explosive  energy  is  used 
instead) ;  however,  there  is  a  critically  important  condition 
associated  with  nuclear  explosions,  which  Hilliar's  scaling  laws 
do  not  address.  Since  nuclear  explosive  shock  waves  travel  many 
times  further  than  those  generated  from  high  explosives  or  TNT, 
and  since  nuclear  weapons  may  be  used  at  high  altitude  against 
flying  targets,  more  general  scaling  laws  which  include  the  ef¬ 
fects  of  a  variable  atmosphere  were  needed.  One  such  set  of 
scaling  laws  were  provided  by  Sachs  in  1944,  and  the  Sachs 
scaling  laws  are  the  most  widely  used  similarity  transformations 
for  shock  waves  traveling  in  a  nonhomogeneous  atmosphere  (see 
Glasstone  and  Dolan,  1977) . 

Sachs  assumed  that  the  total  energy  in  the  blast  wave  is 
independent  of  the  external  conditions,  and  the  only  effect  of  a 
change  in  external  atmospheric  conditions  is  to  change  the  scale 
of  pressure,  distance  and  time.  His  result  was  a  generalized 
similarity  transformation  in  which  the  shock  pressure  and  density 
as  a  function  of  distance  and  time  for  a  blast  with  one  set  ex¬ 
ternal  conditions  can  be  obtained  from  that  for  a  different  blast 


with  another  set  of  external  conditions: 


"If  the  atmospheric  pressure  and  temperatures  are 


changed  to  P'  =  n  P  and  T'  =  0  T,  then,  according 


to  the  law  of  similitude,  the  new  blast  pressure 


p'  (R',  t')  at  distance  R'  =  xR  and  t' 


i  t  is 


given  by 


*>■  (id  ,  f  )  -  ii i’ ( k,  t ) 


where  x  and  t  are  constants  which  determine 

the  change  in  scale  of  distance  and  time."  (Sachs,  1944) 


The  distance  and  time  scale  factors  x  and  x  were  determined  by 
Sachs  purely  from  dimensional  reasoning.  Since  velocity  has 
dimensions  of  length  divided  by  time,  the  velocity  in  the  primed 
system  must  be  related  to  the  velocity  in  the  unprimed  system  by 


u' not,*) 


( H-3) 


and  since  density  is  proportional  to  p/T  through  the  ideal  gas 
law,  the  two  densities  are  related  by 


r  7  i r  r  (*.*■) 


( I 1-4 ) 


When  these  are  substituted  into  the  continuity  equation 
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(H-5) 


and  the  equation  of  motion 


v> 
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the  corresponding  equations  in  the  primed  system  are  produced: 


JO.  If  +  gL  X  Jiv  f  U '  z  O 

TT  Jt:  fr 


( 1 1-7 ) 


3*-  ^  4-  JL  P  ~  o 

x  v*  T  it  i 


( I 1-8) 


?/\ 


where  div"  and  grad"  imply  dif fer entation  with  respect  to  R". 


Equations  (II-7)  and  (II-8)  must  express  solutions  of  the  hydro¬ 
dynamics  equations  (II-5)  and  (II-6).  Equation  (II-7)  already 
has  the  required  form  of  ( I I— 5 ) ,  but  (II-8)  has  the  form  of 
( X I— 6 )  only  if 

X  -  X©  (II-9) 

v 

This  is  one  of  the  relations  between  x  and  x  ;  the  other  was 
derived  by  Sachs  from  a  dimensional  analysis  of  the  energy  equa¬ 


tion.  It  was  assumed  that  the  total  energy  in  the  blast  wave  was 
the  total  energy  of  the  air  behind  the  shock  front  minus  the 
energy  of  the  same  volume  of  air  at  atmospheric  conditions. 
Integrating  each  term  in  the  energy  equation  (p  +-j pu2  +  p J" C^.  dT) 
over  the  volume  behind  the  shock  front  indicates  the  following 
dimensional  relationship  between  the  energies  in  the  primed  and 
unprimed  system: 


9 


E7*')  a  7T  XZ 
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Now,  invoking  the  assumption  that  the  total  energy  is  independent 
of  the  external  conditions,  if  the  mass  (or  equivalently  yield) 
is  increased  by  a  factor  m  from  the  unprimed  to  the  primed  sys¬ 
tem,  it  follows  that 

7T  X3 

Thus  the  two  scale  factors  x  and  t  are  given  by 

i. 
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The  Sachs  scaling  factors,  and  indeed  the  entire  Sachs  simi¬ 
larity  transformation,  can  be  obtained  in  a  slightly  different 
manner  by  transforming  both  the  dependent  and  independent  varia¬ 
bles  in  the  hydrodynamics  equations  (see  Bridgman,  1980)  .  The 
transformed  dependent  variables  density,  velocity,  energy  and 
pressure  represent  dimensionless  ratios  of  the  value  of  each 
variable  to  some  reference  value,  and  the  new  independent  varia¬ 
bles  are  determined  so  that  they  are  scaled  so  as  to  eliminate 
the  dimensioned  constants  in  each  equation.  This  process  of  nor¬ 
malization  of  the  differential  equations  leads  to  the  same 
scaling  laws  as  Sachs',  but  it  has  the  advantage  of  explicity 


demonstrating  that  a  solution  to  the  resultant  normalized  equa¬ 
tions  can  be  used  for  any  set  of  constants  (initial  and  boundary 
conditions)  . 


The  Sachs  scaling  laws  may  be  used  to  produce  values  for 
peak  pressure  (static,  dynamic  or  total)  ,  range,  and  time  of 
arrival  for  a  blast  of  yield  Y  in  a  homogeneous  atmosphere  of 
ambient  pressure  P  and  ambient  temperature  if  the  values  are 
known  for  some  "reference"  blast  of  yield  Y*  and  ambient  condi¬ 
tions  P*  and  T*  as  follows: 


(11-14) 


(11-16) 


It  can  be  easily  shown  (Bridgman,  1980)  that  the  density  and 
velocity  are  given  by 


where  cq  is  the  speed  of  sound  in  the  undisturbed  air.  The 
"reference"  yield  and  atmosphere  is  usually  that  of  a  one  kiloton 
burst  in  a  standard  (STP)  atmosphere  (see  Glasstone  and  Dolan, 
1977)  . 

Sachs  scaling  predicts  the  solution  of  one  explosion  in  an 

infinite  homogeneous  atmosphere,  given  a  solution  for  the  shock 
variables  of  another  explosion  in  the  same  or  another  infinite 
homogeneous  atmosphere.  It  gives  no  information,  however,  in  the 
case  of  a  blast  wave  moving  through  a  real,  spatially  varying 
atmosphere  which  must  be  considered  in  the  case  of  nuclear  burst 
at  one  altitude  when  the  target  is  at  a  substantially  different 
altitude.  To  account  for  the  variable  atmospheric  effects  on  the 
shock  values,  the  nuclear  weapons  effects  community  has  adopted 
the  following  procedure,  called  modified  Sachs  scaling:  instead 
of  using  the  ambient  pressure,  temperature  and  density  at  the 
burst  altitude  in  the  Sachs  scaling  laws,  the  corresponding 
values  at  the  target  altitude  are  used.  Modified  Sachs  scaling 
has  no  rigorous,  theoretical  basis,  but  it  is  very  simple  to 
implement  and  can  readily  produce  results  in  lieu  of  a  step-by- 
step  analysis  of  a  blast  wave  propagating  through  a  real,  varying 
two  dimensional  atmosphere. 


Ledsham-Pike  Alpha  Correction 


Several  investigators  have  studied  the  nonhomogeneous  atmos¬ 
pheric  effects  analytically  (Bach,  Kuhl,  and  Oppenheim,  1975; 
Laumbach  and  Probstein,  1969;  Sachdev,  1972;  Hayes  1968a  and 
1968b) .  Although  these  studies  differ  in  the  assumptions  made, 
most  assume  an  exponential  density  variation  in  the  atmopshere 
and  all  neglect  ambient  pressure.  They  define  various  nondimen- 
sional  parameters  consistent  with  their  formulation  of  the  hydro¬ 
dynamics  equations  and  they  reduce  the  partial  differential  equa¬ 
tions  to  ordinary  differential  equations,  which  are  then  solved 
numerically.  None  of  these  studies  has  been  used  by  the  Air 
Force  for  nuclear  survivability  analysis  because  of  the  restric¬ 
tive  assumptions  made  and  because  they  still  rely  on  an  approxi¬ 
mate  numerical  solution.  An  earlier  study,  however,  has  been 
used  in  survivability/  vulnerability  analysis  (see  Sharp  and  Das- 
sow,  1970) .  This  study  proposed  a  simple  correction  to  the  over¬ 
pressure  values  calculated  or  measured  for  homogeneous  atmospher¬ 
ic  conditions  to  approximately  account  for  homogeneous  effects. 

Ledsham  and  Pike  (1951)  assumed  that  all  fluid  motion  at  the 
shock  front  was  purely  radial  and  that  each  part  of  the  shock 
front  behaved  as  though  the  atmosphere  were  spherically  symmetric 
with  respect  to  the  radial  direction  concerned.  Following  the 
procedures  of  Brinkley  and  Kirkwood  (1947) ,  they  nondimensional- 
ized  the  Lagragian  form  of  the  hydrodynamics  equations,  and  made 
several  assumptions  based  on  experimental  evidence.  Ledsham  and 


Pikers  primary  conclusion  was  that  overpressure  Ap  at  a  loca¬ 
tion  where  the  ambient  pressure  is  P  can  be  approximated  if  the 

d 

overpressure  Ap^  is  known  at  the  same  range,  but  coaltitude 
with  the  burst  by 


AP  - 


(11-19) 


where  is  the  ambient  pressure  at  the  burst  altitude  and  a  is 

a 

the  correction  factor.  The  a  factor  was  calculated  by  Ledsham 
and  Pike  as  a  function  of  distance  by  simultaneously  solving 
their  normalized  differential  equations  for  10,  20  and  40  kilo- 
tons  of  TNT  exploded  on  the  ground.  (Thus  the  a  values  were 
calculated  only  for  the  ascending  portion  of  the  shock  wave) .  The 
values  were  recomputed  for  a  one  kiloton  burst  (no  ground  reflec¬ 
tion)  in  a  standard  atmosphere  by  Alexander.  These  values  are 
shown  graphically  in  figure  (II-l),  and  were  taken  from  Sharp  and 
Dassow  (1970). 

The  modified  Sachs  and  Ledsham-Pike  correction  methods  pro¬ 
duce  approximate  values  of  the  shock  variables  (only  overpressure 
for  Ledsham-Pike) ,  and  for  their  reliable  use  as  predictive  tools 
they  must  be  checked  for  accuracy.  Lutzky  and  Lehto  (1968)  per¬ 
formed  one-dimensional  finite  difference  calculations  of  nuclear 
shock  waves  traveling  in  a  spherically  symmetric  atmosphere  with 
exponentially  varying  pressure  and  density  and  with  an  ideal  gas. 
They  considered  downward  propagating  shocks  only  and  found  a  max¬ 
imum  error  in  modified  Sachs  overpressure  of  about  20%.  The 
remainder  of  this  dissertation  is  a  second  evaluation  of  the 
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modified  Sachs  and  Ledsham-Pike  techniques,  which  uses  a  fully 
two-dimensional  finite  difference  calculation.  The  2-D  solution 
has  the  advantage  that  both  the  upward  and  downward  directions 
can  be  checked,  as  well  as  all  directions  in  between  and  in  addi¬ 
tion,  a  spherically  symmetric  atmosphere  need  not  be  assumed  sin¬ 
ce  a  real,  horizontally  stratified  atmosphere  can  be  modeled. 
Furthermore,  the  effects  of  a  rising  fireball  are  included  in¬ 


herently  in  the  solution.  A  variable  y  ideal  gas  ( 


is  the 


ratio  of  specific  heats  and  is  a  function  of  density  and  internal 
energy)  and  a  variable  gravitational  acceleration  (a  function  of 
altitude)  are  used.  Energy  loss  from  the  fireball  due  to  thermal 
radiation  is  also  included. 


The  particular  set  of  burst  conditions  chosen  for  comparison 
with  modified  Sachs  and  Ledsham-Pike  scaling,  was  a  one  megaton 
burst  at  15  kilometers  altitude.  (This  choice  of  burst  condi¬ 
tions  ensures  that  the  shock  wave  does  not  reflect  from  the 
ground  before  low  overpressure  is  reached,  yet  it  exercises  a 
considerable  amount  of  atmosphere  before  low  overpressure.  We 
consider  low  overpressure  as  15  kpascals.  The  1  Mt  results  are 
not  proposed  as  being  superior  to  1  kt  scaling  -  rather  the 
results  of  scaling  to  large  distances  are  being  compared  to  ap¬ 
proximate  finite  difference  results.)  Prior  to  the  2-D  calcula¬ 
tions  however,  a  series  of  one-dimensional  calculations  were  per¬ 
formed.  To  check  the  reliability  of  the  finite  difference  equa¬ 
tions,  a  steady  state  flow  through  a  shock  and  an  unsteady  flow 
in  a  shock  tube  were  simulated.  Then,  a  1-D  calculation  of  a  one 
kiloton  burst  in  a  standard  atmosphere  was  performed  followed  by 
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a  1-D  one  megaton  burst  at  1  k  kilometers  (uniform  atmosphere)  . 
This  final  1-D  calculation  was  used  as  the  "seed"  for  the  modi¬ 
fied  Sachs  and  Ledsham-Pike  scaling  and  these  results  were  com¬ 
pared  to  the  2-D  calculation.  The  1-D  one  kiloton  and  one  mega¬ 
ton  calculations,  as  well  as  the  2-D  calculations  are  discussed 
in  the  next  chapter. 


h.* 


V. 


17 


i 


**./*./ 


Chapter  III.  Numerical  Calculations 


Finite  Differences  and  Flux-Corrected  Transport 


This  chapter  presents  the  results  of  one-and  two-  dimen¬ 
sional  finite  difference  calculations  of  shock  waves  resulting 
from  nuclear  explosions  in  the  atmosphere.  The  numerical  method 
of  finite  differences  is  used  to  solve  the  hyperbolic  system  of 
partial  differential  equations  of  hydrodynamics  (1-1)  in  one- 
dimensional  spherical  form  and  in  two-dimensional  cylindrical 
form.  The  technique  of  flux-corrected  transport,  and  in  particu¬ 
lar  the  SHASTA  (an  acronym  for  SHarp  And  Smooth  Transport  Algor¬ 
ithm)  finite  difference  scheme  is  used. 


The  origin  of  finite  difference  methods  for  numerically  sol¬ 
ving  partial  differential  equations  can  be  traced  to  a  famous 
paper  by  Courant,  Friedrichs,  and  Lewy  (1928) .  The  method  has 
since  been  put  to  use  in  many  fields  of  applied  science,  and  one 
of  the  most  successful  applications  has  been  in  the  field  of 
fluid  dynamics.  In  this  study,  we  are  interested  in  fluid  flows 
that  are  characterized  by  internal  discontinuities,  or  shock 
waves.  One  possible  solution  technique  is  "shock  fitting",  in 
which  we  apply  special  boundary  conditions  at  the  shock  front 
provided  by  the  Rankine-Hugoniot  equations  (Zel'dovich  and  Raiz- 
er,  1966)  and  then  solve  for  the  smooth  flow  behind  the  shock 


using  one  of  the  excellent  fourth  order  difference  schemes 
available  (see,  for  example,  Roberts  and  Weiss,  1966).  However, 
since  the  shock  surface  is  itself  moving  through  the  fluid,  and 
its  motion  is  not  known  in  advance  but  is  governed  by  the  dif¬ 
ferential  equations,  serious  difficulties  arise  immediately. 

A  slightly  different  approach  was  proposed  by  von  Neumann 
and  Richtmyer  (1950)  which  results  in  an  approximate  method  of 
fluid  dynamics  calculations  and  which  takes  care  of  shocks  au¬ 
tomatically.  Rather  than  incorporate  the  Ra  nk  in  e-Hugo  niot  jump 
conditions  into  the  calculations,  von  Neumann  and  Richtmyer  in¬ 
troduced  an  artificial  dissipation  into  the  equations.  Like 
viscosity  or  heat  conduction,  this  dissipation  spatially  smoothes 
the  shock,  so  that  the  discontinuity  is  replaced  by  a  thin  layer 
where  the  flow  quantities  (density,  velocity  and  pressure)  change 
rapidly,  but  continuously.  This  artificial  dissipative  term  is 
called  the  von  Neumann-  Richtmyer  pseudo-viscous  pressure,  and 
allows  the  differential  equations  to  be  solved  over  the  entire 
region  of  interest  -  including  the  shock  front. 

The  effect  of  pseudo -viscosity  on  the  calculated  fluid  pro¬ 
perties  is  dissipation,  and  by  definition  this  causes  the  proper¬ 
ties  (density,  momentum,  energy  and  pressure)  to  be  attenuated 
with  subsequent  time  steps.  To  prevent  the  dissipation  and  at¬ 
tenuation  from  controlling  the  results  or  unacceptably  influen¬ 
cing  the  results,  especially  after  thousands  of  time  steps,  we 
require  strict  control  of  the  artificial  viscosity  term.  It  can 
be  shown  that  any  difference  scheme  that  attenuates  any  Fourier 


component  of  the  solution  cun  hr  t.  hough »  o  l  nr  .  •  r  •  t  1  i  n  i  rq  v  i  .•.'•ius 
terms,  even  though  no  r.udi  term::,  an.  exp  I  i  c  i  t  1  i-  i  ittor  into  the 
equations.  This  concept  was  exploited  1  v  h  :x  and  Wendroff 
(1960),  and  has  been  widely  used  to  calculate  fluid  flows  con¬ 
taining  shock  waves.  A  major  difficulty  with  the  Lax-Wendrof f  , 
and  other  difference  schemes  using  psuedo- viscos  i  ty ,  is  not  only 
are  they  dissipative,  which  they  must  be  to  allow  the  discon¬ 
tinuity  to  be  smeared  over  a  finite  region,  but  they  are  also 
dispersive  (i.e.,  different  Fourier  components  of  the  solution 
travel  with  different  speeds).  A  distinguishinq,  and  unwanted, 
feature  of  such  schemes  is  the  appearance  of  non-physical  oscil¬ 
lations  in  the  solution  which  unpr edictablv  influence  the  peak 
shock  values  (see  Richtmyer  and  Morton,  1967) . 

In  an  effort  to  eliminate  or  at  least  minimize  these  oscil¬ 
lations,  Boris  and  Book  (1973)  developed  a  new  class  of  finite 
difference  techniques  called  flux-corrected  transport.  Like  the 
pseudo-viscosity  techniques,  FCT  treats  all  grid  points  identi¬ 
cally;  however,  there  are  no  added  artificial  viscosity  terms  and 
no  special  knowledge  about  the  solution  is  required.  An  FCT 
scheme  handles  steep  gradient  problems  by  invoking  a  physical 
property  of  the  continuity  equations  which  Boris  and  Book  call 
"positivity" s 

"An  operator  p  is  said  to  be  positive  if  p  >C, 

j 

all  j,  implies  >C  also.”  (Book,  Boris  and  Hain,  1975) 


In  this  definition 


is  an  element  of  the  vector  (1-2)  and  C  is 


a  constant.  The  thing  that  separates  FCT  from  other  schemes  is 


the  application  of  a  flux  limiter  to  the  antidif fusive  fluxes, 
such  that  the  antidiffusion  operation  is  positive.  These  new 
ideas  are  discussed  more  fully  in  Appendix  A. 

An  FCT  algorithm  contains  three  individual  operations  which 
when  combined  carry  the  fluid  properties  (density,  momentum  den¬ 
sity,  and  energy  density)  through  successive  time  steps,  and  when 
these  finite  difference  operations  are  combined  they  are  both 
conservative  and  positive.  The  transport  operation  advects  the 
fluid  through  the  mesh  (assuming  an  Eulerian  mesh)  and  is  the 
same  finite  difference  operation  as  conventional  schemes  without 
artificial  viscosity  terms  or  other  special  terms  included  to 
handle  shocks.  The  diffusion  operation  is  analogous  to  the 
viscosity  terms  in  conventional  schemes;  however,  the  magnitude 
of  the  dissipation  in  FCT  is  purposely  much  larger  than  in 
schemes  with  pseudo-viscosity.  By  adding  a  large  amount  of  dif¬ 
fusion  in  regions  near  sharp  gradients,  the  dispersively  genera¬ 
ted  ripples  are  overwhelmed.  (Some  FCT  schemes,  such  as  SHASTA 
described  later,  combine  the  transport  and  diffusion  operations 
into  one)  .  The  final  operation  in  an  FCT  algorithm  is  antidif¬ 
fusion.  In  this  step,  equal  and  opposite  fluxes  are  applied  to 
exactly  cancel  the  added  diffusion  almost  everywhere.  The  an- 
tidiffusive  fluxes  are  limited  term-by-term  so  as  to  maintain 
positivity,  and  the  effect  is  that  small  amounts  of  diffusion  are 
retained  near  sharp  gradients  to  eliminate  the  oscillations. 


SHASTA  is  a  one-dimensional,  explicit,  Eulerian  FCT  algor¬ 
ithm  written  by  Boris  and  Book  primarily  from  a  geometrical 
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interpretation  of  the  coninuity  equation.  It  can  be  generalized 
to  include  source  terms  and  thus  allow  solution  of  all  three 
equations  (1-1)  with  the  same  algorithm.  With  minor  adjustments 
of  the  flux  terms,  it  can  be  used  for  Cartesian,  cylindrical  or 
spherical  meshes  in  one  dimension  (with  angular  symmetry) .  By 
alternating  between  1-D  subroutines  in  a  time-step  splitting  pro¬ 
cedure,  SHASTA  can  be  used  for  multidimensional  calculations.  A 
more  complete  description,  including  a  derivation  of  the  SHASTA 
finite  difference  equations,  and  an  outline  of  the  time-step 
splitting  procedure,  is  given  in  Appendix  A. 


One-Dimensional  Calculations 


To  establish  a  data  base  for  the  modified  Sachs  and 
Ledsham-Pike  scaling  of  one-dimensional  results  to  compare  to  the 
numerically  calculated  two-dimensional  results,  a  series  of  one¬ 
dimensional  calculations  were  performed  using  a  SHASTA  flux- 
corrected  transport  code.  A  one-dimensional  calculation  of  a  one 
kiloton  nuclear  explosion  in  a  standard  uniform  atmosphere  was 
done,  followed  by  the  one-dimensional  uniform  atmosphere  analogue 
of  the  two-dimensional  test  case  (a  one  megaton  burst  at  15  ki¬ 
lometers  altitude)  .  These  one-dimensional  results  were  then 
scaled,  using  Sachs  scaling,  to  a  common  reference  for  a  test  of 
consistency.  The  common  reference  used  was  the  one  kiloton  burst 
in  a  standard  atmosphere,  and  this  allowed  the  one-dimensional 


results  to  be  compared  to  the  AFWL  one  kiloton  standard. 

The  first  step  in  any  finite  difference  analysis  is  the  gen¬ 
eration  of  a  grid  of  node  points.  For  the  two  one-dimensional 
calculations,  the  regions  over  which  the  differential  equations 
apply  was  divided  into  discrete  finite  difference  one-dimensional 
spherical  Eulerian  meshes.  These  are  illustrated  in  figures 
III-l  and  III-2  for  the  one  kiloton  and  one  megaton  simulations. 
Since  we  are  primarily  interested  in  the  peak  shock  values,  and 
since  the  peak  values  would  be  constantly  eroded  by  a  coarse  grid 
(large  distances  between  node  points),  the  meshes  were  controlled 
so  that  the  shock  wave  was  continually  maintained  in  the  region 
with  the  smallest  mesh  point  spacing.  When  the  shock  wave  ap¬ 
proached  the  edge  of  the  mesh  (the  criterion  used  was  when  the 
peak  pressure  occurred  within  10  node  points  of  the  last  one)  the 
entire  region  was  conservatively  rezoned,  by  proportionately 
expanding  all  zones  so  that  the  shock  was  replaced  to  the  same 
relative  position  in  the  mesh  as  when  the  calculations  began.  In 
this  manner,  the  shock  remained  in  the  region  of  finest  zone 
sizes  although  all  zones  were  intermitently  expanded.  The  number 
of  times  the  mesh  was  rezoned  depended  upon  the  size  of  the  mesh, 
the  length  of  the  time  steps,  and  the  number  of  time  steps  taken. 

The  initial  conditions  at  the  shock  front  and  behind  were 
taken  from  the  AFWL  one  kiloton  standard  (Needham,  et  al,  1975)  , 
where  for  the  one  megaton  calculation,  these  values  were  scaled 
using  Sachs  scaling.  Using  the  AFWL  one  kiloton  standard  elim¬ 
inates  any  device-dependent  considerations  such  as  energy 


partitioning.  The  initial  pressure,  density,  and  velocity  pro¬ 
files  for  the  one  kiloton  yield  are  shown  in  figures  III-3  - 
III-5  and  in  figures  III-6  -  III-8  for  the  one  megaton  yield.  In 
the  unshocked  region,  ahead  of  the  shock  wave,  the  initial  densi¬ 
ties  and  pressures  were  taken  from  the  1976  U.S.  Standard  Atmos¬ 
phere,  and  the  velocities  were  zero.  The  boundary  condition  at 
the  origin  was  that  of  symmetry: 
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and  no  boundary  condition  (other  than  the  external  ambient  condi¬ 
tions)  was  needed  at  the  other  boundary  since  there  are  no  velo¬ 
cities  or  source  terms  there  (the  finite  difference  equations  can 
be  solved  in  such  a  cell  with  no  change  in  density,  momentum  or 
energy. ) 


The  Doan  and  Nickel  (1963)  equation  of  state  of  air  was  used 
to  calculate. he  ratio  of  specific  heats.  This  method  is  a 
"semi -physical  fit"  generated  from  tabulated  values,  and  handles 
densities  from  10  to  100  times  normal  density  and  temperatures 
from  0.0025  to  1.5  electron  volts. 

The  rate  of  radiation  energy  loss  from  the  fireball  was  cal¬ 
culated  using  a  model  developed  by  McNamara,  Jordano,  and  Lewis 
(1977).  While  the  detailed  radiation  transport  process  may  not 
be  important  at  late  times  after  bomb  initiation,  the  total 
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Figure  III-l.  One  dimensional  1  kt  mesh.  The  burst  point  is 
the  origin  of  the  spherically  symmetric  mesh.  The  mesh  contains 
total  of  100  node  points,  but  the  spacing  is  too  small  for  res^ 
lution  at  the  area  indicated  by  the  bracket. 


Figure  III-3.  Initial  pressure  profile  for  the  1-D  1  kt  calcula 


S*3 


0.00  10.00  £0.00  30.01 

RADIUS, 


RADIUS 


2.5 


RRDIUS.  n  *10 


Initial  density  profile  for  the  1-D  1  Mt  caicula- 


RADIUS,  n  *10 


r  \ 


thermal  power  loss  is  important.  The  volumetric  distribution  of 
this  energy  loss,  as  given  by  the  McNamara  model,  depends  on  the 
total  fireball  power  loss  and  the  spatial  distribution  of  inter¬ 
nal  energy.  The  total  loss  rate  was  taken  from  a  normalized 
power-time  curve  given  by  Figure  7.84  in  Glasstone  and  Dolan 
(1977).  For  a  total  energy  loss  rate  of  P  from  the  fireball,  the 
McNamara  model  predicts  a  rate  of  radiation  energy  loss  from  each 
cell  of 


P 


mi  ui_ 
M  u* 


where  ntj  and  u^  are  the  cell's  mass  and  internal  energy,  and  M 
and  U  are  the  respective  mass  and  total  internal  energy  of  the 
fireball. 

Details  of  the  finite  difference  meshes  and  calculational 
procedures  are  given  in  table  III-l.  The  calculated  pressures, 
densities  and  velocities  are  shown  graphically  in  Appendix  C, 
figures  C-l  -  C-12  at  two  times  for  each  of  the  one-dimensional 
cases . 
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Two-Dimensional  Calculations 
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The  two-dimensional  calculations  were  performed  with  a 
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Speed  (CPU  user/node  pt/time  step)  1100 


time-step  splitting  algorithm  using  a  one-dimensional  Cartesian 
SHASTA  subroutine  for  the  axial  direction  and  a  one-dimensional 
cylindrical  SHASTA  subroutine  for  the  radial  direction  (see  Boris 
and  Book  and  Hain,  1975).  Because  of  large  computer  storage  and 
time  requirements  for  a  two-dimensional  calculation,  the  CRAY-1 
mainframe  was  used.  The  number  of  node  points  was  kept  small  so 
that  the  program  would  run  in-core  and  the  total  central  proces¬ 
sor  time  was  minimized. 


The  two-dimensional  solution  was  made  in  r-z  cylindrical 
coordinates  with  the  assumption  of  angular  symmetry.  Since  the 
motion  of  the  shock  wave  in  the  radial  direction  co-altitude  with 
burst  is  primarily  radial,  and  since  the  motion  along  the  z-axis 
must  be  axial,  the  finest  mesh  point  spacing  was  used  in  these 
directions.  This  causes  the  aspect  ratio  (the  height  of  a  cell 
divided  by  its  width)  to  be  very  high  along  the  right-hand  boun¬ 
dary,  and  the  aspect  ratio  to  be  very  low  along  the  upper  and 
lower  boundaries.  Along  +45  from  the  horizontal,  the  aspect 
ratio  is  one.  The  two-dimensional  r-z  finite  difference  mesh  is 
shown  in  figure  III-9.  As  in  the  one-dimensional  calculations, 
when  the  shock  front  approached  the  edge  of  the  mesh  in  the  radi¬ 
al,  positive  axial,  or  negative  axial  direction,  the  region  was 
rezoned  and  the  mesh  expanded  in  that  direction  (the  three  direc¬ 
tions  were  tested  and  adjusted  independently) . 


The  initital  conditions  were  identical  to  the  one¬ 
dimensional  one  megaton  initial  conditions  rotated  symmetrically 
about  the  burst  point.  The  external  conditions  were  varied  in 
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Figure  III-9.  Two  dimensional  r-z  mesh.  Only  the  upper  half  of 
the  mesh  is  shown,  the  bottom  half  is  symmetric.  The  burst  point 
is  at  the  lower  corner  of  this  figure.  The  areas  indicated  by 
brackets  each  contain  80  equally  spaced  mesh  intervals,  the  spa¬ 
cing  of  which  was  too  small  for  resolution. 
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the  z-direction,  since  pressure,  density  and  temperature  vary 
with  altitude.  The  atmospheric  model  used  to  calculate  ambient 
conditions  is  explained  in  Appendix  B.  The  symmetry  boundary 
condition  (III-l)  was  used  along  the  r=0  axis. 

The  same  equation  of  state  and  radiation  energy  loss  model 
used  in  the  one-dimensional  calculations  were  used  in  two  dimen¬ 
sions  also.  An  additional  source  term  was  included  in  the  momen¬ 
tum  equation  for  the  axial  direction  time  step  to  account  for 
gravity  forces.  A  varying  gravitational  acceleration  was  used  to 
account  the  the  slight  variations  with  altitude  to  maintain  a 
stable  atmosphere  (see  Appendix  B) . 


TJ 


The  two-dimensional  mesh  sizes  and  details  of  the  two- 
dimensional  calculations  are  shown  in  table  III-l  along  with 
their  one-dimensional  counterparts.  The  calculated  pressure, 
density  and  velocity  profiles  are  shown  for  two  times  in  Appendix 
C,  figures  C-13  -  C-18  for  the  radial  co-altitude  direction,  and 
in  figures  C-19  -  C-24  for  the  axial  (r=0)  direction. 


Programming  Notes 


The  one-dimensional  calculations  were  done  on  the  Aeronauti- 
*.V'  cal  Systems  Division's  CDC  CYBER  750,  and  the  two-dimensional 
calculations  on  the  Air  Force  Weapons  Laboratory's  CRAY-1.  The 
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two-dimensional  code  was  vectorized  to  take  full  advantage  of  the 
CRAY'S  pipeline  processor.  In  addition,  the  two-dimensional 
problem  could  not  have  been  done  on  the  CYBER  in-core  since  this 
machine  is  limited  to  about  230,000  octal  words  of  storage. 


The  time  step,  At,  was  variable  and  was  calculated  as  the 
maximum  possible  to  keep  max  {  )=0.25,  where 
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and  v  .  is  the  sum  of  material  velocity  and  sonic  velocity  c  . 


A  significant  price  was  paid  for  a  variable  ratio  of  speci¬ 
fic  heats.  Approximately  43%  of  the  CP  time  was  spent  calcula¬ 
ting  y. 
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Chapter  IV.  Comparisons  of  Numerical  Results  to 
Scaling  Law  Predictions 


The  results  of  the  modified  Sachs  and  Ledsham-Pike  scaling 
of  the  one-dimensional  solution  are  compared  to  the  two- 
dimensional  finite  difference  calculations  in  this  chapter.  The 
three  most  important  parameters  -  shock  wave  position  (or 
equivalently,  time  of  arrival  of  the  shock  wave) ,  dynamic  pres¬ 
sure,  and  overpressure  -  are  used  as  the  basis  of  the  comparison. 
The  errors  cited  for  each  of  the  comparisons  were  calculated  with 
the  two-dimensional  finite  difference  solution  used  as  the  stan¬ 
dard.  For  example,  the  error  in  the  modified  Sachs  prediction  of 
overpressure  would  be 


where  2D  and  MS  subscripts  represent  the  finite  difference  and 
modified  Sachs  overpressures. 


Range-Time  Comparisons 


The  calculated  one-dimensional  shock  wave  positions  are 
shown  as  functions  of  scaled  time  in  figure  IV-1.  These  results 
have  been  normalized  to  a  common  time  and  radius  scale  using  the 
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Figure  IV-1.  One  dimensional  shock  envelopes.  The  shock  front 
positions  as  functions  of  time  are  compared  to  the  AFWL  1  kt 
sta ndard. 


Sachs  scaling  procedure  described  in  Chapter  II.  The  common 
scale  is  that  for  a  one  kiloton  burst  in  a  standard  atmosphere, 
which  is  used  so  that  the  results  can  be  compared  to  the  AFWL  lkt 
standard,  also  shown  in  the  figure.  There  is  a  negligible  dif¬ 
ference  between  the  three  shock  envelopes,  and  the  three  curves 
in  figure  IV-1  are  hardly  discernible.  The  one-dimensional 
range-time  data  for  the  one  megaton  calculation  were  used  for  the 
modified  Sachs  calculations  which  follow. 

The  accuracy  of  the  modified  Sachs  predictions  of  shock 
range  vs  time  is  evaluated  by  comparing  their  results  to  the 
two-dimensional  finite  difference  calculations.  The  slant  range 
R  of  the  shock  wave  at  time  t  after  burst,  according  to  modified 
Sachs  scaling  is  obtained  by  using  the  Sachs  scaling  rela¬ 
tionships  for  a  uniform  atmosphere  with  ambient  conditions  at  an 
altitude  Rcos  6  above  (or  below)  the  burst  point,  where  0  is 
the  angle  that  the  slant  range  makes  with  the  vertical  (see 
figure  IV-2)  .  A  slight  complication  arises  from  the  fact  that 
time  must  also  be  scaled  to  the  new  altitude.  For  the  one  mega¬ 
ton  burst  at  15  kilometers  the  following  equation  was  used  to 
calculate  the  modified  Sachs  range  vs  time  curve  shown  in  figure 
IV- 3 : 

Ra  (l 't)j 

where  P*  and  P  are  the  ambient  pressures  at  the  burst  (15 
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Figure  IV-2.  Geometry  o 
cal  shock  position  is  s 
angle  with  the  vertic 
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the  o  ne-d  imens  io  na  ]  one  megaton  shock  :anq-  calculated  at  It 
kilometers  at  time  t*.  Now  the  complication  a*  ’  -es  because  this 

range  R  occurs  at  time  t  which  also  must  be  coal'd.  Since  we  are 

,  .  * 

interested  in  comparing  ranges  at  a  given  time  t  from  burst  we 
can  take  the  one  dimensional  data  from  those  c:-  -eluted  at  time 


( IV- 2 ) 


where  T  and  T  are  the  ambient  temperature  at  15  km  and  at  (Rcos 

ci  a. 

=15  km)  respectively,  and  we  have  assumed  that  the  quantity  , 

-lJC^is  the  same  at  both  altitudes. 

The  comparison  of  shock  range  as  a  function  of  time  as  cal¬ 
culated  by  modified  Sachs  and  by  the  two-dimensional  finite  dif¬ 
ference  solution  is  shown  in  figure  IV-3  for  three  times  after 
burst.  The  comparison  is  better  for  altitudes  beneath  the  burst 
point  than  for  above  the  burst,  where  the  modified  Sachs  envelo¬ 
pes  are  consistently  higher  than  those  calculated  by  the  two- 
dimensional  code.  The  errors  were  averaged  for  all  shock  posi¬ 
tions  between  1.5  and  5.5  kilometers  for  ^cations  above  and 
below  the  burst  and  these  average  values  are  shown  in  figure 
IV-4 .  Below  the  burst,  the  average  error  in  modified  Sachs  shock 
positions  stayed  between  +5%;  while  above  the  burst,  the  error 
ranged  from  -4%  to  -7%.  A  difference  in  the  scaling  accuracy  for 
altitudes  above  and  below  the  burst  altitude  was  also  observed 
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for  the  dynamic  and  the  static  pressures,  and  the  practice  of 
differentiating  betweeen  the  two  directions  will  continue  in  the 
accuracy  evaluations  which  follow. 


Dynamic  Pressure  -  Range  Comparisons 


If  figure  IV-3  were  used  as  the  base  of  a  three-dimensional 
figure  with  the  third  dimension  being  peak  dynamic  pressure  (or 
overpressure) ,  the  peak  pressures  occurring  throughout  the  region 


could  be  observed  (regardless  of  the  time  of  arrival  of  the  shock 
wave)  .  It  is  possible  to  use  such  a  construction  to  compare  the 
modified  Sachs  pressures  (modified  Sachs  and  Ledsham-Pike  for 


overpressure)  to  the  two-dimensional  finite  difference  calcula¬ 


tions.  This  can  be  done  by  passing  planes  through  the  three- 
dimensional  figure  parallel  to  the  r-z  base  at  selected  pressure 
levels  and  displaying  the  resulting  intersections  of  the  planes 
and  the  figures  on  a  two  dimensional  graph.  These  two  dimen¬ 
sional  plots  are  called  contours. 


The  one  dimensional  dynamic  pressures  are  shown  as  functions 
of  one  kiloton  range  in  figure  IV-5  with  the  AFWL  lkt  standard. 
The  scaled  one  megaton  dynamic  pressures  are  all  lower  than  the 
AFWL  lkt  standard  for  each  range,  and  the  maximum  difference  is 
about  10%.  The  calculated  one-dimensional  one  megaton  results 
are  used  in  the  modified  Sachs  similarity  transformation  below. 
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The  dynamic  pressure  contours  calculated  with  modified  Sachs 
scaling  of  the  one-dimensional  one  megaton  pressures  are  compared 
to  the  two-dimensional  finite  difference  results  in  figure  IV-6 
for  dynamic  pressure  levels  of  50,  20,  and  5  kpascals.  The  modi¬ 
fied  Sachs  contours  were  obtained  by  transforming  the  one  dimen¬ 
sional  ranges  and  dynamic  pressures  to  the  altitude  R  cos  0  , 
where  R  is  the  slant  range  and  0  the  angle  between  the  slant 
range  and  the  vertical  (as  shown  in  figure  IV-2)  .  Equation  IV-1 
was  used  for  the  range,  and  for  the  dynamic  pressure  q,  the  Sachs 
similarity  transformation  is  the  same  as  for  static  pressure: 


( IV-3 ) 


where  q*  is  the  dynamic  pressure  at  R*  in  an  infinite,  homo¬ 
geneous  atmosphere.  An  iterative  procedure  was  required  to  pro¬ 
duce  the  modified  Sachs  contours,  since  the  correct  combination 
of  dynamic  pressure,  range  and  ambient  pressure  must  be  found. 


The  modified  Sachs  dynamic  pressures  are  consistently  larger 
than  the  finite  difference  pressures  at  the  same  range,  and  the 
two  pressures  are  closer  in  magnitude,  for  all  ranges,  above  the 
burst  altitude  than  below.  The  average  dynamic  pressure  errors 
versus  range  are  shown  in  figure  IV-7.  Below  the  burst,  the 
error  was  -25%  to  -35%;  while  above  the  burst,  the  error  was  -8% 
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Figure  IV-6.  Modified  Sachs  and  2-D  finite  difference  dynamic 
pressure  contours  for  three  pressure  levels.  y 
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Figure  IV- 7 .  Average  error  of  modified  Sachs  dynamic  pressure  vs 
range  for  ascending  and  descending  directions. 


Over  pressure- Range  Comparisons 


The  one-dimensional  overpressures  calculated  for  the  one 
kiloton  and  one  megaton  hursts  are  plotted  in  figure  IV-8  with 
the  AFWL  lkt  standard,  where  as  before  the  results  have  been 
scaled  to  the  lkt  standard. 


The  maximum  difference  between  the  scaled  one  megaton  range 
and  the  one  kiloton  standard  was  approximately  12%  and  occurred 
at  a  scaled  overpressure  of  about  700  kpascals.  In  the  compar¬ 
isons  that  follow,  the  one-dimensional  one  megaton  overpressures 
were  used  for  the  modified  Sachs  and  Ledsham-Pike  scaling. 


The  overpressure  contours  calculated  using  modified  Sachs 
scaling  and  Ledsham-Pike  alpha  correction  are  compared  to  those 
calculated  with  the  two-dimensional  finite  difference  code  in 
figure  IV-9  for  overpressure  levels  of  60,  30,  and  15  kpascals. 
The  slant  range  R  to  overpressure  aP  is  obtained  with  modified 
Sachs  scaling  by  using  the  Sachs  scaling  relationships  for  a  uni¬ 
form  atmosphere  with  ambient  conditions  at  RcosQ  above  (or  below) 
the  burst  point.  The  range  R  is  scaled  as  in  equation  (IV-1)  , 
but  now  R*  is  the  one-dimensional  one  megaton  range  to  overpres¬ 
sure  Ap*,  where 


VERPRESSURES  VS 


4P‘  -  4P  A 


( IV-4 ) 


For  the  Ledsham-Pike  method,  the  overpressure  at  slant  range  R  is 


4P  ~  AP 


(■fer 


( IV-5 ) 


where  Ap*  is  the  one-dimensional  one  megaton  overpressure  coalti¬ 
tude  with  the  burst  at  range  R.  The  a  used  was  the  same  as  that 
given  in  figure  II-l.  An  iterative  procedure  was  used  to  compute 
the  modified  Sachs  and  Ledsham-Pike  curves  of  figure  IV-9  to  find 
the  correct  combinations  of  overpressures,  ranges  and  ambient 
pressures. 


Two  observations  are  immediately  apparent  from  figure  IV-9. 
The  first  is  that  there  is  practically  no  difference  between  the 
modified  Sachs  and  the  Ledsham-Pike  overpressure  contours.  The 
difference  increases  slightly  with  decreasing  overpressure  levels 
(longer  distances  from  the  burst);  however,  for  the  levels  shown 
in  figure  IV-9,  it  is  reasonable  to  neglect  the  difference 
between  overpressures  calculated  with  modified  Sachs  and 
Ledsham-Pike.  The  second  observation  is  that  the  scaling  methods 
predicted  overpressures  much  closer  to  the  two-dimensional  finite 
difference  solution  for  altitudes  above  the  burst  than  for  alti¬ 
tudes  below  the  burst.  The  average  errors  are  shown  in  figure 


52 


-  naoiriED  snow 

- LE03HWH»IKC 

.  2-B  3Hft3Tn 


15  kpascals 


IV-10,  as  functions  of  range.  The  modified  Sachs  overpressures 
were  consistently  greater  than  those  calculated  by  the  two- 
dimensional  code  for  a  given  range.  Below  the  burst  altitude, 
the  errors  ranged  from  -20%  to  -30%,  while  above  the  burst  the 
errors  were  from  -6%  to  -20%.  For  modified  Sachs,  this  differen¬ 
ce  cannot  be  changed;  however,  with  the  Ledsham-Pike  correction, 
the  situation  can  be  improved. 

The  difference  between  the  finite  difference  calculated 
overpressures  and  the  Ledsham-Pike  scaled  overpressures  can  be 
reduced  by  adjusting  the  variable  a  .  Using  the  calculated 
values  for  Ap  in  equation  (11-20),  we  can  compute  a  as  a  func¬ 
tion  of  range  as 


Mr)  =  -2 


lo%  A  P  -  / of  A  P 


/«£  Pa.  -  P* 


( IV-6 ) 


This  relationship  was  used  to  calculate  new  values  of  a  for  all 
node  points  beneath  the  burst  for  the  last  36  time  steps  saved  on 
data  tapes.  Using  the  method  of  least  squares,  these  data  were 
then  smoothed  with  a  third  order  polynomial  spline,  and  the 
resulting  a  (r)  is  shown  in  figure  IV-11  with  the  abscissa  scaled 
one  kiloton  standard  range.  Also  shown  in  this  figure,  is  a  plot 
of  the  original  a(r),  which  when  used  for  the  altitudes  above 
the  burst  point  in  conjunction  with  the  new  values  for  lower 
altitudes  produced  the  overpressure  contours  shown  in  figure 
IV-12.  The  differences  between  the  calculated  and  scaled  results 
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Figure  IV-10.  Average  error  of  modified  Sachs  overpressure  vs 
range  for  ascending  and  descending  directions. 


REPRODUCED 


in  the  descending  direction  have  been  reduced  to  approximately 
the  same  values  as  for  the  ascending  direction. 


The  validity  of  this  proposed  new  formulation  of  the  c 
correction  factor  depends  upon  correctly  interpreting  where  the 
errors  occur.  In  particular,  does  the  true  overpressure  contour 
lie  near  the  outside  of  the  "bulge"  directly  beneath  the  burst 
point  (which  would  indicate  that  modified  Sachs  and  Ledsham-Pike 
are  quite  accurate,  as  shown  in  figure  IV-9) ,  or  is  the  true  con¬ 
tour  closer  to  the  value  indicated  by  the  predominant  number  of 
data  points  further  away  from  the  axes  and  the  bulges?  The  avera¬ 
ging  of  all  data  points  in  the  reformulation  above  implicitly 
assumes  the  latter.  The  answer  can  be  found  by  running  the  two- 
dimensional  code  a  second  time  with  a  homogeneous  atmosphere. 
This  two-dimensional  calculation  would  in  effect  simulate  a  one¬ 
dimensional  problem,  for  which  we  know  an  exact  solution.  In  a 
uniform  atmosphere,  the  overpressure  contour  would  be  exactly 
circular  at  a  radius  calculated  from  the  one-dimensional  results 
in  chapter  III  with  the  scaling  procedures  described  in  chapter 


This  second  calculation  was  performed,  and  the  resulting  15 
kpascal  overpressure  contour  is  shown  in  figure  IV-13.  The  solid 
line  is  the  "known"  solution,  and  is  a  semicircle  centered  about 
the  burst  point  with  radius  equal  to  the  range  to  15  kpascal 
overpressure  calculated  for  the  one-dimensional  one  megaton  case. 
Again,  the  bulges  appear  along  the  axes  and  along  the  z=15  km 
line.  When  figure  IV-13  and  figure  IV-12  are  compared,  it  is 


Figure  IV-13.  Fifteen  Kpascul  Overpressure;  Contour  for 
Homogeneous  Atmosphere 


clear  that  the  two-dimensional  solution  is  nearer  the  known  solu¬ 


tion  away  from  these  bulges,  and  therefore  the  correction  factor 
also  is  formulated  with  respect  to  the  two-dimensional  solution 
away  from  the  bulges. 

It  is  reiterated  that  the  uniform  atmosphere  solution  shown 
in  figure  IV-13  uses  the  calculated  one-dimensional  result  from 
chapter  III.  The  3.5  km  range  to  15  kpascal  overpressure  is 
slightly  lower  than  the  AFWL  1  kt  standard  scaled  to  1  Mt  at  15 
km  height  of  burst,  which  is  approximately  3.7  km.  In  fact,  if 
the  3.7  km  range  was  used,  it  would  appear  that  the  outsides  of 
the  bulges  are  the  correct  result.  This  is  believed  to  be  coin¬ 
cidental.  To  eliminate  the  common  errors  between  the  one¬ 
dimensional  and  two-dimensional  results  (especially  the  "clip¬ 
ping"  phenomenon  inherent  in  SHASTA,  which  affects  peak  values, 
see  Boris  and  Book,  1973),  the  correct  one-dimensional  reference 
must  be  the  same  calculation.  It  is  also  noted  that  the 
currently  accepted  ranges  for  the  1  kt  standard  to  overpressures 
less  than  30  kpascals  are  being  revised  in  the  direction  of  the 
results  presented  here.* 


*  Private  communication  with  Lt  Tom  Lutton,  Air  Force  Weapons 
Laboratory,  January  1983. 
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Chapter  V.  Conclusions 


The  accuracy  of  the  modified  Sachs  scaling  method  of 
transforming  one-dimensional  blast  wave  results  to  account  for 
atmospheric  variations  with  altitude  has  been  evaluated.  The 
shock  wave  position,  peak  dynamic  pressure,  and  peak  overpressure 
as  functions  of  position  and  time  as  calculated  by  the  modified 
Sachs  method  were  compared  to  the  results  of  a  two-  dimensional 
finite  difference  calculation  for  a  one  megaton  burst  at  an  alti¬ 
tude  of  15  kilometers.  It  was  found  that  the  shock  wave  position 
predicted  by  modified  Sachs  was  quite  close  to  that  calculated  by 
the  two-dimensional  code  -  the  magnitude  of  the  average  error  was 
less  than  8%.  The  dynamic  and  the  static  pressures  predicted  by 
modified  Sachs  were  both  closer  to  the  two-dimensional  calcula¬ 
tion  for  altitudes  above  the  burst  point  than  for  altitudes  below 
the  burst.  The  modified  Sachs  dynamic  pressure  and  overpressure 
average  errors  were  both  about  12%  above  the  burst,  while  below 
the  burst  they  averaged  about  30%  and  25%  respectively.  Thus, 
modified  Sachs  pressures  for  the  descending  shock  front  (which 
moves  into  regions  of  increasing  density  and  pressure)  contained 
twice  the  error  of  those  for  the  ascending  shock  front.  These 
overpressure  errors  are  only  slightly  higher  than  those  calcula¬ 
ted  by  Lutzky  and  Lehto  for  downward  propagating  shocks. 

The  Ledsham-Pike  correction  of  scaled  overpressures  was 
found  to  be  practically  identical  to  modified  Sachs.  New  values 


of  the  correction  factor  for  the  descending  shock  front  were  cal¬ 
culated  from  the  two-dimensional  finite  difference  overpressures. 
These  new  values  are  the  best  fit  to  the  data  calculated  by  one 


method  for  one  set  of  burst  conditions,  and  to  test  their  accura¬ 
cy  and  reliability  they  should  be  used  in  further  comparisons 
with  other  yields  at  different  altitudes. 

The  SHASTA  flux-corrected  transport  method  was  used  for  the 
one-  dimensional  calculations  in  spherical  coordinates,  and  for 
the  two-dimensional  calculations  with  a  time  step  splitting  al¬ 
gorithm  which  alternated  between  Cartesian  and  cylindrical  coor¬ 
dinates.  In  all  cases,  the  method  kept  the  shock  sharp,  with  a 
thickness  three  or  four  mesh  intervals  wide,  and  produced  no 
observable  oscillations  in  the  solutions. 
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Appendix  A:  The  SHASTA  Finite  Difference  Scheme 


All  flux-corrected  transport  schemes  consist  of  two  steps: 
the  transport,  or  convective  stage  and  the  anti-diffusive,  or 
corrective  stage.  In  the  transport  stage  the  flow  quantities 
(density,  momentum,  and  energy)  are  convected  one  time  step  using 
some  conservative  difference  scheme  which  includes  a  large  dif¬ 
fusive  component.  This  erroneous  diffusion  is  removed  in  the 
antidif fusive  stage  where  special  conditions  are  applied  near 
shocks  to  retain  sharp  gradients. 


The  SHASTA  finite  difference  scheme  used  here  is  the  ori¬ 
ginal  flux-corrected  transport  method  of  Boris  and  Book  (1973)  . 
Diffusion  is  included  inherently,  and  there  is  no  need  to  add 
pseudo-viscosity  as  in  Lax-Wendrof f ,  leapfrog,  donor  cell,  or 
other  "conventional"  schemes.  To  illustrate  the  transport  stage 
of  SHASTA,  consider  the  one-dimensional  mass  transfer  example 
shown  in  figure  A-l.  The  density  profile  is  shown  as  a  solid 
line  at  the  time  Pat  the  node  points  j-1,  j,  and  j+1.  The  den¬ 


sity  is  assumed  to  be  constant  over  the  mesh  interval  (j-—  j+~  ) 

*  2 


The  arrows  at  the  node  points  indicate  displacements  over  the 


time  interval  At  =  tn+1-tn 


The  transported  density  values  f.,  , 


and  fj+i  are  determined  by  conserving  mass  during  the  time 


step  in  a  Lagrangian  sense.  For  example,  conservation  of  mass 


for  the  cell  (j,j+*^)  requires 


(A-l  ) 
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where  S*  is  the  source  term.  (Note  that  in  the  equations  of 
hydrodynamics,  there  is  no  mass  source  term,  but  there  are  sour¬ 
ces  of  momentum  and  energy.)  The  other  cells  are  treated  similar¬ 
ly,  resulting  in 
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To  ensure  that  no  grid  point  convects  past  a  cell  boundary, 
require 


we 


i  Arv 


(A-6) 


This  transport  process  exactly  conserves  mass  in  each  mesh 
interval,  and  if  the  new  density  values  at  time  t**'  were  compu¬ 
ted  by  summing  the  contributions  of  each  convected  cell,  there 
would  be  no  diffusion.  To  introduce  a  diffusive  component,  Boris 
and  Book  linearly  interpolate  between  adjacent  node  points  for 
the  density  values.  These  fluid  element  trapezoids  are  shown  in 
figure  A-2.  Clearly,  the  area  under  each  trapezoid  equals  the 
area  under  the  corresponding  pair  of  rectangles  in  A-l,  and  hence 
mass  is  fully  conserved;  however,  some  mass  that  was  convected 
into  cell  (j--£,  j+£)  is  transfered  to  adjacent  cells  by  the 
linear  interpolation.  This  is  the  origin  of  mass  diffusion.  The 
diffused  mass  is  represented  in  figure  A-3  by  the  shaded  areas. 


Figure  A-3  shows  the  details  of  the  transport  and  diffusion 

of  mass  for  cell  (j-^  ,  j+£).  The  density  for  this  cell  at  time 

M-l  n  TD 

t  is  } •  =  f-  ,  where  superscript  TD  stands  for  transported  and 
0  J 

diffused.  Pj*  is  computed  by  summing  the  contributions  to  the 


cell  from  each  adjacent  trapezoid: 
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Figure  A-2.  Fluid  element  trapezoids.  The  area  under  each  tra¬ 
pezoid  equals  the  area  under  each  corresponding  rectangle  of 
figure  A-l ,  but  there  is  a  net  loss  of  mass  in  the  cell  (j-  ,j+ 


Figure  A-3.  Transported  and  diffused  density.  The  density  at  time 
t  is  the  sum  of  the  contributions  from  the  two  surrounding 
trapezoids.  The  diffused  mass  is  represented  by  the  shaded  area. 
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Now  introduce  two  very  convenient  parameters 


+■  i  ~  (fii  -  <v 


and  substitute  (A-2)-(A-5)  into  (A-7) : 


(A-10) 


Equation  (A-10)  is  a  generalization  of  the  SHASTA  difference 
equation  given  by  Boris  and  Book  (197  3^  .  The  scheme  consists  of 
an  essentially  second  order  differencing  of  the  convection  term 
plus  a  strong  diffusion.  For  a  uniformly  vanishing  velocity 


field,  this  strong  diffusion  amounts  to  operating  on  the  initial 
density  profile  in  the  following  way: 
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For  uniform  zero  velocities,^  =g  ;  but  for  non-zero  velocities,^ 
is-g  plus  small  velocity-  and  wave  number-dependent  terms.  The 
effect  of  these  small  corrections  to  y  are  shown  by  Boris  and 
Book  (1973)  and  Boris,  Book  and  Hain  (1975)  to  be  small,  and  they 
are  neglected  in  this  study.  Thus,^  =  ^’. 


We  wish  to  write  (A- 10)  in  "flux"  form, 


(A-12) 


p  ~T 

where  f  and  f  are  the  diffusive  and  transportive  fluxes,  and 
where  we  have  temporarily  neglected  the  source  terms.  We  know 
that 


-  M  ($-&) 


( A-l 3 ) 


because  the  diffused  solution 


(A-l 4) 
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Now  (A-10)  and  (A-12)  can  be  equated,  and  after  some  algebra  we 
find 


=  [  i tef}1- »,*-■)  *i  ]  f'j  -  [f 


I  ( A-16 ) 


&  ,A~17: 


Written  in  flux  form,  equation  (A-12)  explicitly  states  that  the 

H+l 

transported  and  diffused  density  at  time  t  equals  the  old  value 
at  tn  plus  the  influx  due  to  transport  and  diffusion  minus  the 
outflux.  Appropriate  scale  factors  are  included  to  account  for  a 
non-uniform  grid. 


Another  set  of  scale  factors  must  be  added  to  allow  use  of 
(A-12)  for  a  curvilinear  coordinate  system.  The  generalized 
one-dimensional  radial  continuity  equation  can  be  written  (again 
neglecting  source  terms) 


2l  4-  X  JL  Lvvfv)  =.  o 


( A-18) 


where  V  =0,  1,  or  2  for  a  one-dimensional  Cartesian,  cylindrical. 


or  spherical  coordinate  system  respectively.  Using  f  ^±1  as  the 
in-  and  out-fluxes,  the  finite  difference  approximation  of  (A-18) 
is 


The  one-dimensional  difference  equation  (A-12)  v;ritten  for  varia¬ 
ble  mesh  spacing,  generalized  coordinates,  and  with  source  terms 
is  thus 


This  difference  equation  advances  the  mass,  momentum,  and  energy 
densities  one  time  step  by  conservative  transportive  convection, 
but  it  includes  a  large  diffusive  component. 

The  second  stage  of  flux-corrected  transport  is  the  antidif- 
fusive  stage.  Ar.tidiffusive  fluxes  are  added  to  the  solution  to 
remove  the  erroneous  diffusion  from  the  transport  stage  in  a  con- 
servative,  but  nonlinear  manner.  The  diffusive  fluxes  f...  would 


be  exactly  cancelled  if  these  values  were  subtracted  from  the 
values  everywhere  they  were  added  to  the  solution;  however,  as 
shown  below,  the  a nt idif fusive  fluxes  are  limited  in  magnitude  so 
as  to  eliminate  dispersively  generated  ripples  in  the  solution. 
Although  the  magnitude  of  this  corrective  diffusion  depends  only 
on  the  local  values  of  f  from  point  to  point,  the  antidiffusion 
operation  is  conservative. 

T& 

Since  f  is  the  result  of  diffusive  transport,  the  antidif- 
fusive  stage  of  flux-corrected  transport  results  in  the  final 

values  f  : 
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where  f  is  the  a ntidif fusive  flux,  and  the  scale  factors  are 
temporarily  omitted.  For  exact  cancellation  of  diffusion 


( A-22) 


Such  a  prescription  is  obviously  conservative,  since  every  dif¬ 
fusive  flux  is  added  once  to  the  solution  and  then  subtracted 
once  somewhere  else.  In  addition  to  this  definition  of  antidif¬ 
fusion,  Boris  and  Book  prescribe  the  following  qualitative  limi¬ 


tation  to  inhibit  generating  nonphysical  oscillations  in  the 


solut  ion : 


"The  antidi f fusive  stage  should  generate  no  new  maxima  or 
minima  in  the  solution,  nor  should  it  accentuate  already 
existing  extrema."  (Boris  and  Book,  1976) 

**  A 

The  values  f  are  limited  term-by-term  so  that  no  antidif- 
fusive  flux  transfer  of  mass  can  push  the  density  value  at  any 
grid  point  beyond  the  density  value  at  neighboring  points.  The 
corrected  fluxes  f  are  given  by 


-  s '  i  >  s‘ Vll } 


(A-23) 


where 


Vi  -  C  -  f\ 
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It  can  be  shown  that  (A-23)  is  the  correct  formula  to  quan¬ 
tify  the  perscribed  qualitative  limitation  on  the  fluxes  by  con- 
sidering  all  possible  permutations  of  the  signs  of  f^,  ,  and  & 

(See  Boris,  Book  and  Hain,  Section  IV.)  In  particular,  a  local 

tt>  -  - 

extrema  is  indicated  for  the  value  /V  i£  ei ther  of  /l^a r e 


T  A 


opposite  in  sign  to  f^+j  .  In  this  case  the  flux  limiter  sets  f^yj 


~  A 


-0 .  On  the  other  hand,  when  f  ,  A.’-l  ,  and  fl-,are  a.11  the  same 

n  *  1  *  i 

sign,  the  flux  ]  imiter  leaves  =  f unless  the  j  *“  point 

would  be  pushed  below  the  (j-1)**  ,  or  if  the  (i  +  1)1^  point  would 


be  pushed  above  or  below  the  (j+2)’1-*  .  In  this  case,  (j+j.  is 

corrected  to  avoid  creating  a  new  extrema. 


Using  the  corrected  antidif fusive 
are  computed  as 


fluxes , 


the 


new  values 
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( A-24) 


Time-Step  Splitting 


The  two-dimensional  solution  procedure  was  a  time-step 
splitting  technique,  using  a  separate  Cartesian  and  cylindrical 
subroutine  for  the  z-  and  r-dir ections.  For  a  centered  time 
integration,  provisional  values  are  calculated  for  one-half  time 
steps,  and  these  provisional  values  are  used  to  advance  the 
values  one  full  time  step. 


If  we  represent  the  operator  FfU,v,S,ar,at]  as  that  which 
advances  {u  }  a  time-step  £t  on  a  grid  ir  (or  az  )  using  veloci¬ 
ties  v  and  source  terms  S  ,  then  the  complete  cycle  is  as  fol- 


first  the  r  half  cycle 
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next  the  z  half  cycle 
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Appendix  B:  Stable  Atmospheric  Model 


The  earth's  atmospheric  temperature  profile  was  modeled  as 
consisting  of  six  linear  segments.  Hydrostatic  equlibrium  was 
assumed,  and  the  air  was  treated  as  a  homogeneous  mixture  of  its 
constituent  gases. 


The  exact  values  of  temperature  were  taken  from  the  U.S. 
Standard  Atmosphere  (1976)  for  the  base  of  each  of  the  six  linear 
segments,  and  the  corresponding  pressure  is  calculated  to  satisfy 
hydrostatic  equilibrium 


£  =  -  3  ft)  fa) 


( B — 1 ) 


where  g(z)  is  the  altitude-dependent  gravitational  acceleration. 
The  pressure,  temperature  and  density  of  the  air  are  related  by 
the  ideal  gas  law 


?U)~ 
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where  R  is  the  universal  gas  constant  and  M  is  the  molecular 
weight.  Dividing  (B-l)  by  p  yields 


The  inverse  square 

r*-' 

where  g0  is  the  sea  level  acceleration  of  gravity  (9.80665  m/sec) 
and  r0  is  the  effective  radius  of  the  earth  (6356.766  km). 

Within  each  of  the  linear  segments  the  temperature  is 

TU)Z  \  4.  U  (B-5) 


where  Tk  is  the  temperature  of  the  base  of  the  segment  (i.e.  at 
z=zk)  and  the  temperature  gradient  is  .  Substituting  (B-4)  and 
(B-5)  into  { B— 3 )  gives 

d  (JUT)  =  -  1  - - - - - (B-6) 

H***  >  *  [Tk  4-  U  £l-  *«•)] 

Integrating  over  the  length  of  the  segment, 

A 

where  zfc  and  zt  are  the  altitudes  at  the  base  and  top  of  the  seg¬ 
ment,  and  Pk  and  P*  are  the  respective  pressures.  Carrying  out 
the  integration  gives 
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Using  equation  (B-8),  the  pressure  at  any  altitude  below  86 
km  can  be  calculated.  For  example,  at  the  bases  of  each  of  the 
six  linear  temperature  segments  the  pressures  are  as  shown  in 
table  B-l.  To  calculate  the  pressure  at  any  general  altitude, 
the  temperature  is  found  from  (B-5)  and  then  (B-8)  is  used  with 
(p)  z  instead  of  p  and  z=z  .  With  the  pressure  and  temperature 
known,  the  equation  of  state  is  used  to  calculate  density. 


Table  B-l.  Base  Values  of  Pressure,  Temperature, 
and  Temperature  Gradient 


Altitude,  Pressure,  Temperature,  Gradient, 

km  Pascals  K  K/km 


0 

101325.0 

288.150 

-6.488727 

11 

22702.70 

216.774 

-0.00956044 

20.1 

5446.865 

216.687 

0.9974380 

32.2 

863.9094 

228.756 

2.756184 

47.4 

110.2434 

270.650 

-0.05878049 

51.5 

66.21404 

270.409 

-2.738829 

72.0 

3.841296 

214.263 

-1.956643 

86.0 

0.3737405 

186.870 

0.1087500 

40.0 
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Figure  C-3.  Velocity  profile  fo  1-D  1  kt  burst  at  1  sec  after 
burst . 


Figure  C-4 .  Pressure  profile  of  1-D  1  kt  burst  at  7.67  sec  after 
burst . 


Figure  C-5.  Density  profi]e  of  1-D  1  kt  burst  at  7.67  sec  after 
burst. 
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Figure  C-9.  Velocity  profile  of  1-D  1  Mt  burst  at  1  sec  after 
burst. 


Figure  C-10.  Pressure  profile  of  1-D  1  Mt  burst  at  10  sec  after 
burst. 
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Figure  C-12.  Velocity  profile  of  1-D  1  Mt  burst  at  10  sec  after 
burst. 
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Figure  C-20.  Axial  (r=0)  density  profile  of  2-D  1  Mt  burst  at 
1.168  sec  after  burst. 


>v 


13.00  14.00  15.00  16.00  17.00  16.00  19.00  20.00 

ALTITUDE.  Kn 


\'WV 


ii 


r  ‘Kit  ' 


M  NL-G  N  At  Q  t  3 1  B 1  5  (> )  »  7(  G  {  3  OCi  ,4)  ,  OPUJS  (0:3  0  0  )  ,  SP LUS ( G  S 30 u >  . 

L  h  N  V  r  L  (  C :  3  C  0  )  .  R  (  Q  :  3  b  0  )  ,  7  (  0  J  3  0  0  )  .  -  f  I  (  0  :  5  00  >  •  7N(  0  :  300  ) 

NS  NS  :.CN  :nuT  (3)  ,  JQUT  (3) 

Vi'  S-::C!'  FO  (  C  :  3  0  0  )  .  F  7  (  0  :  3  0  0)  »  7  1 D  (  0  :  3  0  0  )  .  Di/TIJ  (  0  :  3  0  0  )  .  p  3  (  0  :  3  D  P  )  . 

>"  1<  ■:  :  3  0  C )  ,  P  ?.  (  0  s  ?  0  0  )  ,  h  3  (  G  :  3  0  0  )  ,  r-  u  (  G  :  3  0  0  ) 
y  Nr  :c  N  FHO  (  0  :  3  0  0  )  ,  L  (0  :  '5  0  0)  *  GNOME  ( l) :  30  Q)  ,  £  1  (  0!  300).  £11(0:300), 
W  i  &  3  0  C  ) ,  FC  <  0  :  3  0  0  )  .  FON<  0  :3  00  >  ,  E  £T  A  (  0  :  3  00  )  .  FN(  J  !  30G>  ,F£<  0  *  300  ) 

,  >  (  0  :  .3  J  C ) 


i 

-  N  ^  M 

1  j  =  Q  ,  7  n: y  - 1 
J  )  =  0 . 

?.  J^0iv-l 

L  L  £  <  ..I )  -  C  . 

L  U  L  (  J  )  =  G  . 

?.  :  =  i ,  4 
j  , : )  =.o . 

- 1 


.1N'£  IN  Z.T  I  A  LI.Z£ ..Th.L.-.Mlil.Sti . . . . . _ . .  ...  . . . .  .. 

L  L  I  NT.  AT  a  (  A  ,  N.  M  ,  NM,  I  OUT  ,  JOU  T  ,  NRO  KS  ,  NCOLS  ,  V  ,  TST  OP  ,NSTOP  ,N! 

e T ,  t i m e  i n :tl p., r , z ) . . * . . .  .  . . 

U ’ CrN5 ' £F 

'  1-v  :  MT:  a  L..  .C0NCIL1.0N.E . . . .  . . . 


!.  L  Cl."r  [.  AT  A  (  A  ,  N  ,  N  .  NM,  lOUT  ,  JOUT  ,N^OWS  ,  MC  0  L  £  ,  D  T  ,  T  I  M  £ ,  NS  T  £  P ,  P. ,  7 . 0 ) 


T  ./  2. 


AN'.:,;  th;  PCOPEPTIcS  x  N  TH  l  FACIAL  CI'-lCTION. 

DO  P.C  J=0,y-1 

■;T  -  h  ■;.  F  -HALF ....  S.J  ,iLP. .  . . . .  . . . . 

N'C  N.'NC':'Y,  f-  A  NO  Z  “MO  M  2  NT  IJ  M  ,  AM!,  c  Nrl  CGY  ONE -HA  IF  STEP. 

O  AL  L  7Cl.CE  M'  (  a  ,  N  ,  NM  ,  J  ,  CEN‘7£L,  HT  ,u-  ) 

C  AL  !.  OLE  TIP  (  N  .n-T  ,  Q°LUS  ,CC  N7EL  ,  F  ) 

C  A L  t  L  GE  T.C  tr.L  i .S.PL US].. . . . . . . . . .  , . . . . 

call  C  ha  ST  AR  (  NM  ,7  (0,1)  ,  0  ,  A  ,  N,  J  ,  CPI  US  ,  S  PL  U  S ,  R  ,  FO,  F  T  ,  7  T  0 , 3  7  T  0 
p  r  .  L  ,  P  2  ,  P  3  ,  R  A  ) 

CAL.  L  CsiACAR  (  N.M.,7  [  C,3  )  ,2  ,  A.,  N,  J  ,  OnL  US  ,  2  PL  US  ,  R  ,FC,  FT  .VTC.DVTl 
17  F  ,  ~  l  ,  i:  2  ,  p  3  ,  A  t, ) 
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CALL  SHfi£T/!P(KM.V(C.2»  *1.  A.  N.J,  CALL'S*  S  PLUS*®  tFCtPT,VTD.DVT,Jf 
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CALL  S  EE  ’  m  Z<  A  ,  N  ,M. ,  NM.t  I. » H.T  i  S.PL  UE  *  Z)  .  .  .  .  . . 

r  all  Shaft  AZ  (KP.  V<0.3  »  .2*  A,N*  M,  :,QphJS.SPLU3  ,Z,PO,FT,VTO, 
c  rVTr,FC,Pl.(:2> 

CALL  SS£T£ZCA*h»1*M.I.-7  iCPlUL.Z) 
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r 

C  AL  L  V  "  l  Fc  ?(  A  »N»Ht  N3,  ;  ,  V,  -  *0,  r  ,  r.wnu*'  1 ,  t  n  *  WS  ,  F  0 .  PQN .  T  A  , 
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CALL  V  E  L  Fc  Z ( A  »N*M»NM»I»V/.  FHO.  F ,  GMON£.£  1 .  £  1 1,  WS  ,F  0  ,  FON  ,  BET  A  . 
C  FN,F£,EX) 


si  is  7  itut:;  nil  values  into  a. 

CALL  SU2Z(P,\/.hJlL*N.'l.l.)-... 
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lFfNFAXF.GE.N-10)  CALL  RE  ZCNF.  (  0  .  A  ,N  ,  m  ,  N  M  ,C  ,  Z ,  FN  ,  7  N ,  V  ) 
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IF  (NPAXZ^.GE  .fi-lU)  CALL  REZO  NC  ( *  1  .A  .N.I.Ni.  R.7.  RN.ZN.V) 
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C  CJlCULJT  j  ciCI/'LLY  CENTERED  VELOCITIES  rOP  POW  J 
C  _ _  _  ....  . . 

CO  10  IC=l.N-2 
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10  CONTINUE 
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^  .  I  3= 5*NH+ J*N  + 1 

fr.  I4=5,NM*J*R 
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E  .  o?=  (MIC  +  ll*5  CCM/2. 
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13  CONTINUE 
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SU  PS  T  I  T  LT  _  f.CNSITY*  R£  DIAL  A  ML  AXIAL  MOMENTA.  AMD  f  NEPCiY 
V  SCRATCH  A-r/Y  INTO  A.  FOR  RCW  J. 

r.O  1C  CC-C.N-2 
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a  (  n  c :  r.  >  =  \i  ( :  c ,  1 ) 
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13  c  ont  :nu-: 
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~  Nr 


111 


J  J  -•  J*  ^ 


:  4 ' )  / : c 
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'll 
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II^ENSEON  £(Qt7*NM-ll.\/<0:*-1.41 

SUBSTITUTE  LENSITY.A.XDU _ AMD  -  FACIAL  MQUZNT A.  ANIL OJ £££.!.. F£DM_Iil£- 

SC  PATCH  A-FAY  INTO  A,  FOR  COLUMN  I. 


CO  12  JC  =  C  .M -  2 _ 

N  OOF  =  I  +  JC*N 
A  <  N  0  C  t£ )  =  V*  C  JC.i) 
A(NH4N0CC)  =_V.(J.Cj  2). 

A (2*NM  +  N0l  E)  =  V<  JC.3  ) 
A(3*NM+NCCE>  =  V(JC,4) 

CONTINUE  .  _ _ - 

P ET  UP  N 
E  NC 


SU3F  CUT  INF  i/ELPFP  (A.»N»  MM.  J  .  V».M  .RHQ  ,L.  GMQNE.E  J.£ll»  WS.fQ. FQN. 

C  £ £T  A  »  F  N  ,FftEX) 

CIMENSION  £( Q 1 7*NM-1) , V ( 0 1 M-l , 4) 

S  IMS  NS  ION  RHQ  ( 3  1  N-ll  .£  Ml  N- 11  .  GMOM  ( Q 1 N  -  1 1  .-E-lIQi  Nr  11  .Ell  (2  liLrl.L*  — 
C  WE  (  Cl  N-  1)  ,  FO<  C  IN-1 )  ,FQN(  0?  N-ll  .  °£TA(0  :N-1 1  ,FN  (0  IN-1)  .  FE  (  01N-1  ) 

C  ,  E  X  (  C  *  N -1 » 

LOAC  FA  CAL  VELOCITIES  ANC  PRESSURES  INTO  A  EMM  SCRATCH  AR- AY  V  FOP 
-  Ow  j  . 


J  N  =  J  *  N 

1  1  =  4  *NM+ JN 

2  ?=6*NM*JN.  . . . .  . . . . .  .  _ 

;  ^=  £,  *NM+  JN 

co  ic  ;  =  o.  n-  i 

1G  toil)  =  VC  *2)  iC  VI  1*1).. . -  . .  . — 

ro  ic  :=c,n-i 
15  A  (!♦  ]?l  =  V  (If  )  )/Vd  *  D 

co  2o.  :  =  o,.n- i  _ _ _  .  _  - -  -  - 

23  A  c  +  ll  )=CVHOT  (C  .  .A  CC  H  ♦fiaS(A(Min.LM.:-l> 

00  25  1=  0, N-l 

A  O  I?)  =CVfGT  »flC±JL.2).*A?SlA!.I+I?)  ?  «J.T  ,1.,  E -11  ..  _ . . . 

V<I.2)=CVNCT(0.,'/C»2),A)S(VC,21).LT.l,E-3» 

25  VC,3)=CVHCT(0.,V(:,3),APS(V(I,3)).LT.l.E-3> 

CO  3  C  I=0,N-1 _  ..  .  •  - - 

E(  I  >=A>3S(<VC:,4)-VC,2)»A<:C  11/2. -W<I.3)*A(I*I2  1/2. )/»/(:,  1)1 

30  RHCC)  =  V(I,D 

CALL  AIC  IF  NO  *  t  *  1  j  2.Si  3  M0  N.£ » 1 »  .  E  F  »  n  ,  i  ,  £  1 1  ,WJI,  FO.F  ON, QLT A  ,FN  ,.Ff  *£XJL. 
no  4  C  I  =  Q,N-1 

40  A  ( I  *  1 3  )  =  CMON-  (I  )*=  ( I  )*PHC  (I  ) 

OETUF  N  . .  .  .  .  - 

t  Nr 


SU3R  CUT:  KE  M1LP221A  ,N*.fcL,  NM.«  L.  V/ .RHO  .L.£MQNEE»  E  J.*L1  i.  WS.-FQ* -ECLN, - 

C  3E7A,F!v,"E.EX) 
dimension  ac  q:7*nm-d  , \n  o  im-i  ,4> 

L  :M£  M  E  C.N  f.HQ  (J-tM-riJ  *.£aQXM-l)  ♦O^ONE  taiH- 1)  .,-EXXiUM.-  1)  ♦£±-LU-iM-l  )  « . 

C  WS(oSM-l)  ,FO(QiM-l)  .FONIOS  M-l I , BETA ( 0  I M-l »  ,  FN  (  0  !M-  1)  *FE  I  0  l  M-l ) 
C  ♦EXIDtr'-l) 

C  LO  AC  AXIAL  VELOCITIES  AND  PeESSUFLS  I N”  0  A  FROM  SCRATCH  ARRAY  V  FOF 
C  COLUMN  I. 


z  1=5*NH+: 

: 2=6* nm+: 

I  3=h»NM+.I  ..  _ _ _ _ 

CO  1C  J=Q,*-1 

i z  a  (n* w+:n  =  u j,2)/tf  <j, i > 

CO  15  J=  0,^-1  _ _ _ _ _ .  _ _ 

ir>  MM\t:?i:dij,.ji/vu,ii 
CO  2  C  J=0»F-1 

2U  A  ( N  *  1 1.=  C >/ M G.T.CJJAJ  AJLN_*J.±I  U.t  AQ Z  LA  W.JAl  1) ,U_L  L^JL  _ 

HO  ?5  J  =  0  » 1 

A(N<»j*:2)=CVMC1(0.,A(N*J+I2»,APS(A(N*J*-:2)).LT.1.E-1) 

*  <J,  C»  =C  VK-  l Q *_» V 1  J^. IMU sJl ’  >  •LL*1»E.“3> _ 

25  »7(J,  -)=C\-MGT{C.,V(J,3»  ,ABS(V(  J,3)  ).LT.l.E-3» 

CO  ^(]  J=0,K-1 

c.  H  C  (  J  >_=  V  i  J  .  1  ). _ _ _ _ _ _ 

3  0  f  ( J)  sASS  (  - V<  J»2  >*MN*J  *1 1»  /  2.-V  (  J.  3»  *A  (N*  J  +  12)  /  ?.  »/\M  J.  1>  ) 

CALL  AIC  «FFO,F,  1,293*GM0NE»1.  f  6  »M  ,  -  1 ,  E  1 1 ,  W5  .  F  C  ,r  ON  ,  EET  A  ,  F  N  ,  FF  ♦  FX  » 

CO  AC  J.=  0 % t'-.l _  _ _ _ _ „  _  _ 

40  A  »J*  N*I3  )  =CMONt  J)*RHC  (J  ) 

R ST  URN 

END  _ _ _ _ _  _ _ _ _ 


SuaSCL’lINE  M  A  XT  (A«N*£L*NM»QT  ».l  *  Z  *NMAXS.*NM  AX  ZL  »  NMA  XZh)  .  _ 

Clr'El'CI.Ch  A(GJ7*NM-1)  ,R(  0  IN- 1)  ,7  (CJM-ll 

r:uc  *Ax;nj »  \,'lcc:ty  anc  calculate  time  etep. 

DT=1.E1Q 

NUM  2=  f  M- 1)  /Z*  N .  . .  .  .  ..  . . 

NU*  1  =5*NM  +  MJm2 
NIJM  =  i*NM  +  M.M  2 

PMAX-A(MJM)  ...  .  _  _ _ _ _  _ _ _ _ 

on  ic  :  =  i,n-2 

0^1  =  0,  ?*5  MMI+1>-P<I)  )  /  (  A  6E  (A  <NUM1>:  U  +S  l}P  T  (  1. 4*  A  (NUM+I)  / 

C  A  (NIJM 2 ♦!.).)  ).  .  . . . . . . 

: f  c  r - 1 ,lt,ct) r7=cTi 

:  F  (  A  ( MIN  ♦  I  >«lT  , PM AX ) 00  TO  10 

Pm;  >r  a  (NUP+I  *  _ _ _  .  . . 

N  M  A  X  ~  =  I 
10  CONTINUE 

P*A  X  =A  (NL'M  )  . . .  ....  . . 

NUy  1  -f-HM 
Nt.JM  =  t*NM 

CO  2  0  J=  (M-l  >/2+l.tM-i  ....  .  . . 

OTl=  C  ,  (  7(  ,j  4- 1  >  -7  (  J)  )  /(  ARE  (A  (NUMl  Lj’NII  +EQ5T  (1.  4*  A  (NUM*  J*  N  )  /  A 

C  (  J*N)  )  ) 

:  c  ( c  1 1 »  l  T » c  t  >  c  Tj=.c  j  i _ _ _ _ _ _ _ _ _ 

I  F  f  A  tMJf'*J*N>  ,L7,P*AX)GO  TC  20 
P  MA  X  =  A  (  NIJM  ♦  J*  N  ) 

NMAX2H=J  . .  ..  .  .  ...  _ _  .  . . .  . 

20  CONTINUE 
PMAXrA(MJN) 

(?o  3Q  j=i,  n- 11/2  . .  .  . 

DTI  =  0 , 25  ♦  (  2C  J>-7  ( J-l>)  /  (  A  dZ  (A  (NtJM  HJ’NII 
C  +SC-T  (  i  ,4*  A  (NUM+ J*N  (  J*N  )  )  ) 

iFicTi.LT.mcTsru 

T  PC  A  CNUM  + J*N)  .LT.PMAX  >  no  TO  .10 
P  MA  X  =  A <NU*+J V  > 

NMAX2L-J 
30  (CNTINuj-- 
i?  Et  Lp  N 


LU?r-.CL!'r.M  RJ.20K£  HOIf.  .A.  U«  M,  NM  .  -  .Z.r»N.7N,  y/l  .  -  -  -  - 

l.'!*r  NSICk  l  C  0  8  7*NM-1»  .MfllN-l  )  «  7 ( 0  8M- 1 ) »  kN ( 0  8N-1 )  ,  ZN  (Oil-l) 

c  v  n  :y-i  .*♦> 

912  =  19-1  >/2  . 

F  N(  C  1=0. C 
ZN<Mi?)=ieeoj. 

IFCCM  2  SC.  75 

P^ZONJ  T  p  i.  t-fi  LF  OF  T rl z._  MLS.P  • 

r?  rcNTiNt'E 

CO  Zb  J=  M 1 2+  1 1 W-1  _ _  -  - 

26  ZN(J>=(7(J>-Z(Ml?n*l.li+Z(M12) 

wMf_  W  =  M12 

JL«  5T  =  M12  ...  ..  .  . . 

7NM  =  ? (Ml  2) 

CO  4C  J= 9 12,9-2 

C  o  2  <?_  :  =  C  1.N- 1 _ _ —  . . . . — 

on  29  K=1  ,9 
29  V  f!  ,*<|  =  0.0 

7NP  =.  <  Z.N  1 J)  *ZN  f J  +  l ) )  /  _  „  .  _ 

70  Z01DF=  (7  (JKFk>*7  ( JN£W*1)  ) /?. 

: <m  7np . r . zolcpj  tfcn 

JN'£  W  =  JNE.W+  1  _  _ _ _ _ _ 

:fc jnew.ec  .*-i>  ~hfn 

CALL  £  T  AT  9  8-l»J»A,N,M,N9,RN»ZN,Ml?8 

GO  -r  41  . . .  .  _..  .  _ . . . -  - 

if  NO  IF 
GO  T  C  7  0 

7  2  =  Z'JM 
C-N  =  7NP-7NM 

CO  37  JJ=0»JN£W-Oi.AST _ 

71  =  CM  JiriT  +  JJ)  *7  (JLAST*  JJ+  l)  >/2. 

:  (M  wJ.  rO.  JtF.W-JLASTI  Z1=ZNP 

C  =  Zl-72  . . . .  ..  - - - 

CO  36  I=l,N-2 
CO  76  KK  = 1  ,4 

lb  7  <  I  *  K  K  »  -  V(1  ,KK»_tA(.l.K.K-l  )  *  NM  >1  ♦  (  JUACT*  J  J)*N  )  5.G^DFN _ 

72  =  71 
37  CONTINUE 

JLAST  =  JNtW  . . . .  . _ . 

7N9  =  7NP 
V/  <0*1)=V«1,1I 

\l  (0  ,2)  =  V  (1  ,2  > . .  . .  .  ..  .... 

1/  (0 , 7)  =\l  ( 1  ,3  ) 

\J  <  C  ,4)  =  V  (1  ,4  I 

\M N  -  1  ,1  )  =V.fN'-2.  *1> _  .  .  .  - - - 

V(N-1,  ?)  -\i  (N-2,2) 

V  (N-  U  7  )  =  V  (N-2,3) 

V(N-J,4)=V(N-2,4) 

CO  3  (  I  =  0  , N-  1 
NOr.f  =I+J»N 
UNOr-flrVCttl 
A  (  N  m  ♦  N  O  C  f  )  -\J  (I,  ?) 

A(2*NNUNOCE>  =  V<:,3) 

A  (7»NM4N0CrM  r  y  C  ,  4) 


A  (5*NM+N0C£)  =VIL,21/V . 

A  (E*KM  +  N0CE)=V(T*3)/tf  Ctl) 

NUMtfF =5* hH+NCOE 

NUH  VZ  =  G*NM+NQCi. _  _  .  __  __  ....  _  .. 

E  =  (  V  (I »  4  )  -  <A  (  NUM  VR  1  *A  ( NUN  WEI*  A (NU^tfZl^A (NUMtf  Z) 1/2. »/V l 
CALL  GAS (V C ,1) ,E,1.293.GM1,1 .F 61 
A  t4*  TNfNCCE)  =GNl*E*y  11,11  .  ....  .  _ 

30  CONTINUE 

40  CONTINUE 

41  t  CNTINUE  ...  ..  . .  . .  _  .  ._  .  ....  _ 

JO  4?  J-N12M.M-1 

42  Z (J)-ZN(J) 

r'.CTU-N 

PsZ-TNt  in  tfe  -  ac.  al  cirlctiof 


50  CONTINUE 

CO  51  :=1,N-1 

51  F  N  ( I )  =  F  f  I  1  *1 . 11 . . _ . 

I NE  W  =  1 

:last=i 

F.NP  =  <EN(01  +R.N  tl.)J.Z2a . . .  ..  .  .  _ _ 

CO  55  1  =  1, 5-2 
CO  54  J=  0  »  F-  l 

CO  5  4  K  =  .1*4 _ _ _ _  _  ._  _ 

54  V  <J,k)-Q.  0 

R  NF  =  (RMI  )  *PN  (I  +1 )  1  /2  • 

55  ROLCF.?  r^JI5ikLLt£JIN£»Lt4lI/2.  .  ......  ...  . 

I- ( P5P.0T  ,-OLCP)  THEN 

in-:w=:new^i 

IF(If£rt»£G  »Nr.il_  .T.HtLN _  _ _  .  .  . . 

CALL  S~ATM3,I,A,N,M,NM,RN.Z  ,mi?) 
go  tq  be 

£  NO  IF  _  _ _  . 

GO  tG  55 
ENC  :p 

F  2  =  .F.NK  .  _ _  _ _ 

0eN  =  i-  NF*?NP-PNH*PNM 

m  he  :i=o»:new-:last 

*1  =  IE  (ILAST+II1  HR  (ILA5T  +  ::  +  11  )  /  2 «  . 

IF(II,  Et],:Mrt-:LAST»  R  1  =  -  N E 
F=T1*4  1- f  2  *F  2 
CC  hi  J=  0  I  F-i._ . 

co  ei  kk=i,4 

ei  \l  (J  ,  KK  )  =  V  (J,  KK  )  +A  UKK-  1>  *  N*  +  {  IlA:/  )  ♦  J*N)  *  F/OE  N 

F.2  =  FI  . .  . . 

52  CCNTINUE 

:la:t=:new 

c  nh  =  fnp  _ _  __ 

co  6 :■  j=:,n-i 

nocf=:*j*k 

A  ( N  0  L  E  )  =  V  (  J  t  1) 

A  (NM+NOCE1  =V  ( J.  2) 

A(?*NM  +  NO<~e)-\MJ»3) 

A ( C  L  fc )  =V ( J, 41 
A (5*NN+NOC 5) = V ( J. 21 /V ( J, 1 ) 

A  <fc*NM*NOCn  =V(J,  3)/v/<  J,  1  ) 

MJM\/P  =  5*Mr  ♦NQLi. 


1U9.  NUM\/  2=£*NE  ♦  NOCc _ _  -  _  -  _ 

110.  t=l  vMJ.4  )-  (A  (NUM\/R>  *A  (NUMVRI*  A  l  NU  M  7  )  *A  (  NUMtf  Z  )  >/2.)/VCJ,i) 

111.  CALL  GAS  (\MJ,1>  ,E. 1.293, GM1.1 .561 

112.  A  (<4*NM  +  NCC£)  =  GMl*£*v/l  J,ll_  -  ...  _  . 

113.  63  CONTINUE 

yf.  114.  65  CONTINUE 

115.  EC  CONTINUE  ..  .  ..... 

116.  DO  6  7  1=1,  N-i 

117.  67  t  (I  )  =  RMI) 

1  1  ».  J  E  T  'J  F  K  .  _ _  ...  .  _ . .  .  ... 

r 

f  r.  7  on-  "h>  t  C  V  t'  -  HALF  OF  TH-  HESH. 


119. 

75 

CONTINUE 

1  PC  . 

CO  7  F  J=l.*l? 

131  . 

76 

7N(M12-J)  =Z<*12)-  (2  (N1.2)  -2  (M12-J)  )  *1.11 

122. 

JN-W  =  M12 

123. 

JLAST  =  »-12 

124. 

7.  NN  -1  (M12 1  _  _ _ _ 

125. 

00  90  J=  G  .  E'12-2 

^  £ 

TO  7  <  1=0. N-l 

127. 

QO  7«.  K=l.‘» _ _ _____ 

126. 

79 

V (I ,K)  =  0.0 

1  2C  . 

ZNP  =  (7N  (P12-J)  +ZN<M12-J-1  1)  /.?. 

130. 

6Q 

ZOLCF.  =  (Z  (JNEW)»7(JNEW-lM/2  . _  _ 

131. 

IF(znp.L-'.7ClCP)  teen 

n?. 

JNEW=jnew- i 

1  33. 

I  f  <  J.NE  w  i  £  Q ,  Q  )  _IhiI  N  _  ..  . . .  .  .  . .  . 

1  34. 

CALL  StATM1,J,A,N,M,NM,RN,7N,M12> 

1  35. 

GO  Tf  9 1 

1  3C. 

ENC.IF  _ _ _ _ 

0137. 

GO  fc  c  0 

1 

z  NT  Z  F 

1  3°. 

22  =  ZM*  _ _ _ 

l-.O. 

CEN  =  7NM-ZNF 

1  -41  . 

CO  6  7  JJ  =  0  .JLAl’T-JNEW 

1  -.2. 

Z1  =  (  Z  (  JL  AS  ”  -  JJ)  t2.(  JL.AS  T  -  J  J-  1 )  )  /  2 . 

143. 

IF( vJ.EO. JLAST-JNEW)  Z 1=  ZNP 

1  u4  . 

0  =  72-71 

145. 

co  it  :=i.h-2_ . . . .  . 

luC. 

TO  6  E  KK  =  1  ,4 

147. 

86 

SI  (I  ,KK)=V  (I.KK)  *A  ( (KK-1»*NE+:  ♦  ( JLAST-JJ)< 

14°  . 

Z  2  =  Z 1  . . 

1  4C, 

67 

CON  TJNUE 

1  5C  . 

JLAST=JNLK 

151. 

ZM-  =  ZNP  .  _ _ 

15  2. 

tf<G  ,1>=\M1,1) 

142. 

tf (0  ,  2  »  =  V  Cl ,2  ) 

lr  4. 

sl  (3  ,3)  =.V  (1  li) _  _  .  .  .  . 

1  5C  . 

si  (0  .A)  =v/  (1,4) 

lr  5. 

\MN-l  ,1)  =  SI  ( N -  ?  .1) 

15  7. 

V  (N-1,2)  =  V  IN-2,2)  . 

1 . 

tf  <K-l,3)=V(N-2,3) 

144. 

V  ( N  -  1 , 4)  =  V  (N-2,4) 

1‘  G. 

DO  6  c  :=G,N-l 

1  1  1  • 

'Zi  _ 

NCCT  =  I*(E12-J)*N 

N) *  G/DE  N 


1  3  , 


A  (NCC-)  =  Vd»l> 

A  (N'M.NOCF  )  =  Sl  C  .?) 


1%4.  A  (2*NM*NQC£)  =V  <I„3i 

u5.  a  (3*m»kcce>-v(I«  <•> 

A«6*NM+N0C£)=\M:.2>/\M:*1» 

1 67.  A  t6*KM*N0CE»  =V(I*3)/>/ Ch»1  I 

1(8.  NUNVR=5*K^*NCCt 

If  5.  NUMV7=fe*Nf'+N0CE 

1/C.  E  =  (  V(I,4  )  -  (A  (NUM\Z£)  *A.I  NUMA/f.)t  A  (KUMtfZ)  *A(NUMi/ZJJ/2.)/ VII.1J 

171.  CALL  GAS  (V  (I  « 1 1  «E»1.293*GM1»1  .EE) 

1/2.  M4*l'f'*NCC£)  =  CM1*E*\/(I  » i ) 

173.  B“  CONTINUE  .  ..... _  ._  ...  .  _ _  . _ 

174.  90  CONTINUE 

179.  <^1  CONTINUE 

176.  CO  92  J=1.M2  . . .  . .  .  .  _  - _ 

177.  9?  7(M12-J>  =  7MM12-JI 

1  7 c .  -  F.T yc n 

179,  E  NC  . .  _...  _  _  .  .  _ 


1 

? 


£  UH  r.  CU  !  I  K£  S  T  AT  f*  t  ILIk  ,  IJ  .  A  ,  K.  fi .  M .  \  .  Z  •  Ml  2 ) 
CTMEM^QN  A(  D  :7*NM-1I  ,Tv(  0  J  N-l  )  ,  ?t  o  iM-1) 


IF  l  ILIF.)  2  5, 10.75  .  .  . 

M-SH  THt  UP*»tF  ►*  CUMjfl:  Y  FF01  4  XI A  L  flCD'  U  '0  1-1 

25  CONTINUE 

CO  3  G  J=IJ,M-1 

CALL  ?■ or  l/l  j  )  ,r  ,^,P)  . . . 

CO  3  C  1=0,  T-  1 
NCCE=I*J*N 
A  (N0C£)  =  C 

A  (NM  +  NOC  E ) =0  .  0 
A  <2*NM*NOCH  =  0.  0 
A (3*NM+N0C£)  =  £ 

A <4*FM+N0Ct)  =  F 
A  (5*t'M  +  NCCt)  =  Q.  C 

A  (6*NM+N0Ct)  =  C.C  .  . . 

30  CCNT!M,£ 

5  F  T  U  R  M 


M:  SH  THf  ►  I G  T  jCUNDA-Y  pc01  cAOIAL  NOt.E  IJ  TO  N- 


53  CCNT.IN.Ut  . .  . . 

oo  55  :=: j,n-i 

CO  5  5  J=0»M-1  . 

CALL  PfCF(Z(J),Cj£»P) 
N  or  t  =  I  J  *  N 

A  (NO  CE  )  =  C 

A  <NM  +  NOC  FJ  =0.  0. . 

A <?*NM*NOCf ) =0. C 
A  (3*NM+NOC£)  =  t 

A  (A*NN  +  N0CE)  r_.p. 

A  (5*NM  +  NOC  t»  =0.0 
A  (e*NK  +  KCCE)  =  0.0 
5?  CONTINUE 
ALT  urn 


WISH  T  m£  LCfcEF  -3  0UNDAF.V  F  5.0.1  AXIAL  NOT  I  j  TO  0. 


C  ONI  T  I N  U5 

no  sc  j=:j,mi2 

CALL  F--CF  f?(Ml?-j)  ,0,t  ,P» 
DO  RC  :=C,N-1 
N  OC  £  =  I  -f  ( Ml  2  •  J  )  *  N 
A  (NO  CC  )  =  C 
A  (NM+NOC E ) =0 . 0 

A  (2*NM*NQCE>  =C.D  ..  . . . 

A  (3*N1  +  NOr  £)  =  F. 

A|i«*M«-NCCE>  =  P 

A  ( 5  *  NM  4-  N  DC  L )  =  0.  Q . 

£  1 6  *  N  N  +  N  O  C  * '  >  =  0.0 
C  CNTXNUF 
P  F.  T  U  P  N 

r  N[ 


SUORCL’INi:  P-CP(ZM 

ciMt-rsio  2(o:ic>* 

JATA  7/ 

♦  0  •  •  • 

♦  .i.74  0  30CCCZ*02.  . 

♦  ,94C0u03GC£  +  G2,  . 
DATA  1/ 

♦  .26€1GC00C£+G3,  . 

♦  .2  7C650P00E+  03,  . 

♦  .1677  4  COO  OL +  03.  * 
DATA  P/ 

♦  ,ldl325030£+C7,  . 

♦  .ilC242330£+ Q4,  . 

♦  «°  0  45 59  56  GE  *00*  . 
PO  =  6  3  5  f  ,  ?  6  f 

GO  =  9.8  06  65  . 

t-  =  >*.3143 
HT1  =  ?a.9644 
IK»  =  ZM*  0 . 0  0  1.  „  . 


,j.L  .PASCAL  ) 

T  (0  *  10)  ,P  (0  S  1  C> 

11QQ00G  OCE+  G2. 
515  0C30  00L+  02* 
980  C  00  0  J  0  E+  02, 

2167740  0  0E+  03, 
2704C9GQQE+ 03, 
1917  200  0  QE+  Q  3  , 

227  C  2bc  69E  + Ob, 
6621404Q2E+  03, 
4496376090  JO  , 


. 2 JlflOOQ JOE* 32, 
.72Q0003Q3E+  02, 
. 13Q0CQQQ0E+03/ 

.  2166^7  QOOE  +  03, 
.  21426333QE  +03, 
.195080 JGQE+Q3/ 

. 5446:65395+  35, 
.  364129611E+  02. 
.  319224329E+  00/ 


.  322  QQQQGQE+92 
.  660000  00  0E  +  02 


.  22  6  756  0  Q  3E  +  0  3 
.  186670  0  3  02  +  0"* 


.  86  390  9  Li 5F  +  0  4 
.  373740  «*t  OEL.+  C  1 


CO  1Q_.:=Q,10 _ ..._ _  .. 

IFCZkP.LE.ZO)  GO  TO  11 

10  CONTINUE 

11  CONTINUE  _ _ ..  ... 

SLOPE  =  ITO-IC-DI  /  (ZO-ZC-l)) 

A  =  -GG*R0*R  C*WT>/K/SlOPE  /  (RO +Z<:-1) -T  <1-11  /  SLOPE  ) 

5  =  <Z.CI- 1  )ri KMJ/.ic  Q.+2KH.L/  <F  0.+  Z  (X  - 1  n .  _ 

C  =  Aior  (  (F0 +ZKM) +T  c-ll  /  (C  0+ z  (  :- 1  )  »/ (T  ( 1-1)  +SL3  PF¥(  7KH-7  ( 
C  =  C/(cO+7<;-l)-T<l-i>/SLOPI  ) 

PQYNE  -  P(:-l).*£XPCA*_LdrC) )  ._ 

temp  =  T  C-l) +SLOPEMZKM-Z  (I- 111 
PASCAL  =  PCYNE^O.l 
C=PASCAL,I»TM*Q  »  QJJ 1/R/  T  lUP 


-n  * 


••  L  =  PASCAL/ 1.5/0  .  . 

P  =  4  .  E  r 

GM1P  =  1,5 
NO  =  1 

If  T  =  (EL+t>>/2. 

CALL  GAS  (C,t  IN'-,  1.29  3  ,  GM  1  ,1  ,£b> 

I  F  f  AES  (G  Ml  -Gm  ip  LL*lL*  £-9)  GO  "  C 
I F ( NC. GT . 1 £0  »  GO  TO  50 
2 P  =  PASCAL/n/GPl 
IF(lP.Gt.£:Nt)  ,£L  =  eint 
:f(ef.l',,e:nt)  fh  =  -:int 
NC  =  NC+1 

CHIP  =  GM1  . . . 

on  tc  20 
f  =  f:nt+c 


1. 

2  . 

7 


4 


l 

( 

c 


r 

[■ 

7 

t 

9 


ID  . 

1 1  . 

1?. 

13. 

1  4  . 

c 

r. 

r 


17. 
i «. 
ib. 
2  C  . 
.71. 
22. 
.??. 
74. 
??•-. 
D  . 
77. 
2*. 
7«>. 


'  D  . 
<1  . 
<?. 

<3. 

5  4. 

3  5  . 

V  . 


SuPii  CU'If-E  R  ACI/T  ( A  ,N«.M,  M  H.  E.  7 . 7IH£,.0T)  -  _ _  _  ~ 

CIHENSIC  H  A<  C  :7*NM-1)  .Rea  XN-1  1  ,21  rjSM-l>  ,P<1^  >  *T(  ie» 

C  A  A  P/.04»,279,.(:«.c4,l.,.$7/:',.7,.41..297..192..i42».115».0<*7» 

C  .  C72.  .ce.  .  052,  .032  ,  .022.  .32/  _ _ 

PAT.  T  /  0  .  ,  .25. .5, .75.1. ,1.25,1. 5, 2. .2. 25.2. 5  ,2. 75, 3. ,3. 5, 4. .4.5, 
C  £  .  , 4.  ,  1 C.  / 

CALCULATE  TCTAL  FIFDALL  N 3  f  Y  LOT. 

HOTSLM  =  Q.Q 
CO  2D  J=0.H-1 
CO  2  C  1 =  0  »  N-  1 

NCCE  =  I+J«N  .  _ _ _ 

clNT  =  A  (3*NTNOCL)/A  ( NO'  E )  -  (A  (  5  *  N  M  *  N  0  G  E  )  *  A  (5*  N  M  *N  0  C  E  )  ♦ 

0  A(!:.*N4+NOOE»  *A  <6*  NM+NOCE)  )  /?, 

IF(£:nT.Lc.7.£5)  GC  T0  2 r>  .  ..  _ 

VOL  =  3.  14l5c/>-.*  <  <c<  I +1) +MI  >  )  M  (D1  > +- (I )  >  -  (R  <  I-l  >  ♦R  (  D  >  * 

C  (F  C-l  )+R  c  )  )  )  *  (ZC  J  +  l) -Z  (J -1  > > 

hC’SLh  =  HCTSUf..  ♦  A  (NOC£>  *  VOL  *EINT  *:.INT*EINT  ..  _ _ 

20  CONTINUE 
25  CONTINUE 

CALCULATE  EACH  ELtNE  NT  *S  ENERGY  LOEj. 

THAT  =  (T1MDCI./2..  )  y...Q..a3i  .  .  . 

00  ?C  1=2.19 

I  F  (  T  Cl)  . G b  .THAT)  GC  TO  31 

30  CONTINUE  .  ...  _  .  ...  . 

31  CONTINUE 

I  r  <  I  .LT  ,  1R  )  THEN 

P HA T  =  P  O*  (T(:)-:hA1)*  (P(I) -PC-1)  )/lT  ID- Hl-i)  ) 

c  l<?  c 

FhAT  =  C.C2 
END  IF 

P”0 T  =  PHAT  *  1 « 491 13 
CO  4 5  J=0,H-1 

CO  4  C  1=3,  N-l  .  ....... 

NOPE  =  I  +  J  *N 

EINr  =  A  <  3  *  N  M  ♦  N  O  0  E  I  /.  ( NO  )  -  {A  (r*NN+NOCL  )  *  A  (5*NM+NOOE)  ♦ 

C  _ _  _  a<  <:,*nh+nO0E)  *  A 16*  nm  +  nocE)  )  /  2, _ 

I F( EINT.Lt  .7. E5)  GO  To  45 

O  =  RT  flT  /  H  CT  3  UH#  A  ( NOG  t  )  •  f  I  N  T*  r  ;  M  *  ;  NT*  O  T 

A  U^NP  +  NOCE)  =  A(3*NM  +  NOEE)-q  .  . 

40  CONTINUE 
4r’  CONTINUE 

c  £  *  U c  N  .  ...  . .  . . . 

r  N" 


N*.* 


121 


E  U  5 C U T I  N  £  I  SCiTA  (£*U*!1#NM»  IJU7  •  J0U7  »N:;.nWE*MCtXLS  »  V  «TSIOP*-  -  — 

C.  NS-'OF.NflTr  ,  0  T  .  T  IME  ,  K  S’  f  P  ,  7 ) 

0IM£  hEICN  /  (  2  J  7*NM-1)  *  IOC  T  ( .3)  ,  JOU"  (3 )  ,  'J  (  C  IN  -  1  *  4)  ,  R  (0  SN-1 )  .  Z  (  0  SM-1 
NET  E  F  =  0 

£  A  P  (5,100)  "STCP.NSTOP,  Nt  :T£  ,NFE.  N_T 
100  F0C*1AT  (£10.4,41?) 

READ  (5,200  TIMt.NP.  OWE  »  NCQL  E 

£  £A  C  1 5  *£5  0  G  >  ( IOU  T  ( I )  .  1  =  1 .  N"  ON  E  )  .  (  JOUT  C),I  =  1.NC0LE) 

2  00  rc^ItTl  flQ.^.315) 

3  00  FO:.«/T  (£15  )  .  ..  .  . .  .. 

i  P(NrS  .EQ.  1)  GO  TO  11 
c  F  A  C  ( F  )  UM.NETiP.DT 

FEAD13)  IF  ( I )  *  I  =  0  •  N  -1 ).  •  ( Z  ( J ) »  J  =  0 »  H  - 1 1 ...  _  _ 

0  0  1C  J=C,M-1 

10  P:AD<8>  JC,(A(I*M-*J)  ,i=3,6) 

G  0  T  C  2  0  . . .  . . . . . 

11  CONTINUE 

12  F  E  A  C  ( 8  '  TIM  ,  NET-iP,  o- 

P.EACta)  (R  (I)  »I=Q.«N-1.L»(Z.(JI,  J=C*N-1>.  .  ...  _ 

CO  1 3  J=0,M-1 

13  f  F  A  C  (  8  )  JC»(  MI*NM  +  J>  ,1=0,6) 

IF  (NETIF  .NF.  NSTJ...GQ  TO  12.  .  .  _ _ _ 

20  CONTINUE 

GO  3  C  NOC F  =  C,N-1 

GO  3  C  .1*1,1  ..  .  _ _ _  _ 

3  3  \M  N  0  [F  *  I  )  =  A  (<I-l>*NM+NOCE> 


SUE  •:  CUTINc  OLTCATAIA,  At , M«  NN  ,T  OU  T  ♦  JQUT  ♦  NM  OWE  ,  NCUL  S-»£T  *  . - 

C  TIHF.K'ITFF,*  ,Z,:STOP> 

r  :h  : ks; c n  1 i  a  :7*nm-d  ,  :out  ( o  ,  jout  < .? )  ,  p ( a  :n- d  ,z  ( o  :m-i > 

WE  :  IS  (h,  ICO  TIf*Z.NST£P,u' 

1  06  FO-  P  XT  <//////,  13V  ,  »T:  Ki.=  t  ,  lu  .  4 , 5  < ,  *  NS  *  P= »  ,  1  *3 .5  X  ,  rOT  =  t ,  *.  1  j  .  4 ) 
WEI'1’?.  (N,7Q0>  T 1  K  f  i  N'ST  E  Ft  u* 

700  PC- ,MJ5T  <E  lb. 5  ,I1C.£15. fc ) 
on  i  ;  =  a,rs-i 
l  a  r- 1  T  ;  ( - , o  c )  PC) 

cc  2  J-C.f'-l  ...  .  .  .. 

^  WFIM^aci  7(J) 

CO  j  J:|],NM 
Z  0  7  1  =  0  , 6 . 

1  ^R"t(9,80l!)  Mint'tJ) 

?  0  J  r  0:  Mf.T  C15.E  ) 

CO  4  0  JC  =  1  ,NC  CL.I 
KCT  d  C>, 2  0  0  JOC"  (JO) 

■’03  PCP^fli  (i  nx  ,»CCLUMN=t,  I  5) 

;  it,,  3  0  0..  ..  ..  .  .......  . .  .  .  _ .  _.  _  J 

■SCO  F0F*irrf//.5X,»Jt,6X,tZ(J»t»llX,tn)-0t,6>,t‘R  -i  01EN  T  UN  t ,  3X  ,  +7  -MOMENT] 
CMt,4X,tcNEFGXt,ev,»PR:SSUf<Ct,bX,tF-y/:L0C;TY^,5x,tZ-VEL0::TYt) 

[0  15  JsQ.M.  ....  .....  .J 

ZHT  =  2(J)-Z((H-i)/2) 

15  wPITElb,40C)  J.Z^-T,  (A  C*NH+J0U"  tJO+J*N)  ,1  =  0  .6) 

4  00  F  0?.  1  AT. .( 3  X  *  13  ,.£X3X*£.liUJ»l  J  ■■  ...  ...  . . . 

4  0  CONTINUE  ’ 

CO  70  I0  =  1,NCCWS  ; 

Wr.I  '£(6,5  00  IQJCr.LIQ)...  ; 

ECO  FORMAT  ClCX,*-CWst, 15)  j 

4  p  :  T  £  (  h  ,  f;  0  C  ) 

t  30  FOF.  (//,=>  ,  tjt,  6X.*_f_£..(  JA 1 .» 11 X  ,  f?  nQt.bX,  ?JLrM  OiQLT  UK?  #  3X  » tJL-ir.M£Nj 
O  *  ,  4X,  t£  NE  FG  Y  t,  6X,  t  PP  E.  CSU  R  £♦,  5X  ,  t  -  VELOCITY*,  5X,  *Z-VEL0C  IT  Y*> 

cc  4E  :=  0 ,N-  1 

45  WPIT  S  (b,  43  0  I  tc.  C)  .  (A  (  J*UN  +  It  I0LT  CO)  *N)  .  J=D  .  6)  . . . . 

7  0  C  CN  *  :  M  '£  , 

1  F  (  :  C'T’oF  .  C  C*  i)  )  P  •_ T  U c. N 

wc:t^(7)  t:hl.nu£P.,l.t  .  ..  _  . 

W  c  I  ”  :  (  7  )  (F(I),I  =  Q»N-1),  (Z  (J)  »  J  -  0  ,  M  -  1 ) 
no  <3C  J=0,AV-1 

EO  W  ^  I  T  E  l  7  )  J,(A(C*NM+J)  ,:=0  ,f'l 


7o  ccn~: 
if<:c 
wc:t^ 
w  =  l"  : 


S IJ r.  C LI T  :  N A I  R  (f.  hC .  IX  . RhO Z . GM  ON £ . L  Z ♦ N , 1 1 . £11 . WS*  F 0  .F.ON,  3LI.A  *.  _  ... 

C  F N , rr  ,  _  1 

: r  k  fhc  ta  :n-i>  ,ex(  a  s k-i  i  . z <  o  ;n-i > , gmone  i  os  n-u  , eh  q  :n-u  . 
c  £i  i  (0  :k-d  ,wsi  q:n-d  ,f u  i  o:  n-ij  .ronim  n-u  , git  ai  q*n-.u  .fniuin-ij 

0  ,  f  i  (  Q  S  N  - 1  > 

r£  Ml-Ph  VS  ICAL  FIT  TO  'HE  EQUATION  CF  STATE  Of  AIR  _ _ 

sugfout:n,  ;y  l.<-.  loan  and  g.h.  nickel. 

7  r  M  FJ  £  ~  A  T  U  E  t  S  F£CM_,.Q25  TO  1.5  LL-LtRON.  tfOLTS.  _  _ 

tENSITIES  FROM  ia**l-7)  NORMAL  DENSITY. 

P  £  1 5  U  r  £  *  (GAMMA-1 .)  *F.HO*£  ,  Wh-S  GAMMA  IS  A  FUNCTION  OF 
D£ N SI T Y  /NS  ENERGY. 

F.r  0  =  matl?IAL  ItNSITY. 

Pl 2 * * * * 7-  OZ  =  1.2R3-:t  ^FGAG-AMS/CU  •  IC  CLCM-  Tt  -  .  IN  THE  UNITS  OF  T  ML 

FECGLEM.  .  ...  . .  - 

...  -  -N  £R  GY/ MASS 

i?  =  1  Jr.SK/Mlir.AC-rAM,  IN  THE  UNITS  Pr  T  >•£  PROIiL  EM.. 

C- v  ON  £  =  (GAMMA  -  1.) . 

SO  1  1  =  0.  N-l 

1  ill)  =  EX (I ) .  — . . 

SO  2  I  =  D  »  N-l 

2  r  (It  =  A('S  (E  (  I  )  )/■-:? 

SO  3  1  =  0  ,  N-l  . .  .  _ _ 

3  El  (;>=<"  .c-£  (I  ))/ Q  .975 
CO  5  1=0. N-l 

P£l  =  0.,575*  IFPCII)/'MQZ1**Q.Q5  .  . . . 

I  '!  =  »•.•'♦  Q  .  35  7*  ALOG  (RHO  C  )  /-  HO  Z  )/?.  3  0  26 
5  EllCIMf El-tCD/Cll 

CO  IQ  I=C,N-1  ..  ...  .  ... _ _  _ _ _ 

10  WS  (I >  =  1. / (RXP (-£11 (II  > +1.  > 

ro  2  0  :  =  o ,  n-  i 

20  WSII  )=C(/rGT(0.,WS  (If  ,£1(I>.L£  .-5.  »  . . 

L.0  3  C  I  =  G  •  N-  1 

30  W  S  ( I  )  =  C  V  VG  t  ( i . ,  vt'S  ( I )  ,  £  1  ( I  )  .  GT  .  5  ,  ) 

SO  A0  :=  0  .N-  1  ...... 

R  0< I ) =WE  C)*->P(-r.  C)/A.46» 

u  o  f  cn  ( : )  -  ( i .  -  w  i  ( i ) )  *  -  xp  (  -e  ( : )  /c. .  r  3  > 


1  O  5  c  _  =  Ci .  f-  1 

50  SET/. (I  )  =  (  C.U  i£*WS  (I  )♦  C  .3  3  2*  (1  ,-wi  C)  )  >  *  A  LOG!  A  MAXI  ILII!  .1.  )  >.  _ 

C  /  2  •  3  0  ?  E 

[O  EG  I  =  0  ,  F  -  l 

2  f .  2  =  A  •  *(»l-cc;>  /.F.H  CZ  )  *  *  0 . 0  3  5 

E  £ ■.?  =  <•*>,*  C-hOdl/CHCTI  **0.  0157 
!  2=  <*_ (  .) -£F2>/L"2 

Vi  WS  ( I  )  =  1 .  /  (LXR  {-£2  »  +  i.  1  ._  ...  . .  . . 

so  sc  :  =  c.k- i 

P.5  i  1C  »  =  (  <:  )-■*£.  1*0. 33  3  3.33  3  3  33  *?  ?3  3  3  3  3  3 

CO  7  0  1=  0  »N- 1  ...  -  ... 

7  J  WS  (  I  >  s'-tfMAT  (  j  .  ,vr  ( I)  ,  £  1<  1  )  ,L£  .  -  5.  > 
cn  8  0  I  =  l;  •  N-  1 

wsc  »  =  LVMGT  (  1  hi.  (!)  .  £1(1  )  .GT  .  5.) 

UO  R  0  1=0.  N-l 

F N  ( I )  =  *5  ( I  )  *  I  >P  (  -  E  1 1)  /  25  .  5  ) 

OT  -il’A.  (I)  =  2  ETA  C)  Ml. -WE  (I)  )  t  J.  J45*  vS(I) 


ru  loo  :  =  g»n-i 
t.  3=  (  c  ( :  )  -  It  0  .  >  /  6 . 

10C  r  L  (  :  )  =  C'y  MT,  T  {  t  . ,  1.  /  <£XP  (-£3  >  .LT  .  3.  > 


CO  no  I-C.N-l 

FhC  F£C  =  tf  PCt  I  I  /  FhQ  Z  )  ♦*.!•£  7  A  f  I) 

110  GMO  N;  (H  =  (  C.  It  1+0  .  2‘5£j*FO<  I  >+3 . 2t*r  ON  ( I  )  +  3  . 1  3  7 * P'J  ( I  )  ♦  0 
C  )  *  R  i".  C  F  £  C 
m.EtL!FN 
L  NO 


S! 


SUjRIUTINL  Sr.  All. HZ  (MM  .V,I  ,A,N  .  V ,  I C 0_  .G PL  US.  S  PL  LS  ♦  Z  .FC  .  FT  .  ID  .  .. 

c  c\/tc.fc.ci.k2) 

OI^Et'CICN  A  (Q  J7*NP-1)  *  VM  3  t  M-l  >  ,  CPUJO  (  0  5  M- 1 »  t  EPL  J  3  l  0  S  M-  1  >  .  Z  (  0  » P-1  ) 

C  ,  CC  (0  JP-i>  ,FT  (  OJJi-D^j/LC  (0  SM-1)  10  !M-  1)  ,►  C  I  flZM-1)  .  _  _ 

C  h1( 3t“-l) ,F2»CSM-1) 

i  nh  =  :  *  h* 

ii=:i'K+iccl  . 

: 2=:i+n 

ADVANCE  Thl  FLUIC  F.P_DP.£R.T  ONI  'IM-.  SUP  IPI.  AJCIAI _ DUiCTIQN. 

FQ»  C  C  L  L  y  N  I C  CL  .  P£  TURN  NEW  \/AlUES  IN  V  . 

: F  :  =  0 .  CECITY  . 

1,  FACIAL  ,J,C4_  N'  tH 

2,  AXIAL  P  CM  ENT  UM 

3  .  ;N  .CCY 


m  ( o)  =  a  <: ucoi > 

\i  t  c  <  c )  =  o .  c 

PU  5  J=l,y-L 

c  1 «  J'  =  (7  ( J  )-7  ( J-,1 '  >/.<Z.U*  i  >-7.  (  J-t  >  ) 

F2fJ)=(7(J  +  l*-Z(JM/(7(Jft)-7(j-in 

5  CONTINUE 

f C < 0  >  =  J,  0 

F  T  (  0  )  =  0  ♦  0 

00  1C  J=  1 ,  p-  H  .  . . . . .  .  . . 

Fn<j)-o.?E*(A (n*j  +  :i> -a(n*j+: ?)  \ 

FT  ( J ) = ( C  PL  IS  ( J ) *QPLUS  (J)-2.*OPLL'j(J>«'0.7'5)*A  ( N* J  +  1 1 ) 

C  -  <GrLLE  U>*op.LU£  U)  -Q*2.51*a<p*J*.2)  . . . . 

10  continue 

CO  2  c  1 1  v- c, 

VTr(w»=A(T»J*:i>+Rl(J)*(F0CJ-l)4F-(j-m-R?(j)*(FD(J)+FT(J)) 

C  -  2.  *D?  (o)  *QPLUS  f  J>  *EPLU5  (J  )  1  (J)  *  (1.  -  QPLUJ  (  J-l) ) *EPLUE  (.1-1 

20  ca^INOC 

VU  <  t1'  —  .■(  )  r  \,T  (M-L) 

r  vu  (c>  =  n.  c  ...  _ 

JO  ,3  0  J=l,P-4 
J  vT  C  ( J  )  =  v  *  C  (  j  ♦  1 )  -  v/T  1  (  J  ) 

30  CONTINUE 

orn  (m-7  i  =  c.  o 

F  C  i  0  >  =  U .  0  _  .  ..  . _ 

n  0  40  J=  1  ,  )<-4 

S=rifN(l.,CVTC(J)  ) 

A 1*  A  *1 N1  I S  •0.VTja.Jrl.l/_t.2!  4 »  /  ?.  »*«  F£  <rv_rD  (J  )  ,»*fl  ♦  IZS  ,  ...  _ 

C  S  *C7Tr  ( J*  n/F?  U>/2  .)  j 

fc  (  j  )-  an  ex  i  <  o  , ,  Al ) 

43  CONTINJE  ..  . . . .  .  .  ..  ......  _ 

ro  sc  j=  i»p-<« 

V  <  J>  ( J)*CC<  J-1>*IU.U>  Vt-FC  (j)  »L?(  j)  *2,  _  „ 

OQ  C  CN  T  I  N'Ut  . 

=  A  tl  N  +  IC  C3L-*-<  M— 3  >  *  N 1 

1 1-2)  -  a  c  ^♦icolm  V2>  *ni  * 


v/M-n  =  L  Cl^+ICCLt  IM-l)  *N) 
IFCI  .  £0 .  l.AN3.ICQL»£0»  L>  TH£N 
LO  6  C  J=  C  » f'-  1 
t.G  '/  C  J  )  =  Q  . 

£  NC  IF 
(■  '  (:  T  U  F  N 
F  NT 
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C  LJ n  r  CUMM  Sr-.ASTAF.  (NM,  V,K  ,  A  .N  ,  J  ,  QPLUS-,  CPLUS.f  ,FD  ,  FT .  VT  D.  L  V-4  D 
OFC.Tl,F2,r3.FA) 

niKSNSICN  £<  0  !7*NM-1)  . V( 0  tN-l  )  .  CPLUI  ( 0  S  N -1 )  ,  SPLJS  (  0 *  N-l  I  ,  F  (  0  « N-l, 

0  ,FL  (Q  :N-1  )  ,f  T  <  Q:N-1).  Y  7i_  <3  IN-1)  .IVTCICJN-  11  .FC  l  OJN-1)  . 

C  1  (  0  IN-  1)  .  f?  (  J  IN-1  >  ,R3  (  0  IN  -  1  >  ,  F  u(  0  :  u-  1 ) 

T NM  =  K  *  NM 

JM  =  J  *  N  _ 

1 1=  ;nm+ jn 

c  .  . . .  _  __j 

(.  £0  V  A NC  £  ~!-E  FLUID  FFQPEMY  <  CM:  'rI^-  STEP  IN  Y  HE  FACIAL  OIFECTION; 

C  FOW  J.  cjtupK  tM  N  E  ^  VALUES  IN  Y.  ; 

C  IF  K  =  0,  CLNS1T.Y  IS  ADVANCED  _  _ 

f  1*  FACIAL  MOMENTUM 

r  2,  AXIAL  MOMENTUM 

C  3,  IMF  GY 

C 

r 

SO  5  I  - 1 1  N-M  . . . .  . 

F  1  <  I  I  =  <F  (I  >+F  (1-1)1  /  tt/F  (D 

F2d  *  r  ( f  ci^i:*i)i/2,/sii) 

f-3<  I)  =  (F  (Il  -E  lI.-l))  / _ U  C  +  1>-FC-1))  _  _ 

faC)  =  (c  c  +  D-r  mi  /  (t(i+n*n:*iii 

5  CONTINUE.  ; 

c  ....  . . .  .  _  .  ...  _ _  _ 

F  c  <  0  )  =  0  .  C 
f  M  0  )  =  0  .  E 

so  ia..:=i.N-A _ _ _ _ _  .  .  . . .  .  _ 

FT  (I  )  =  Q.25M  A  (I +11) -A  C  +  I2)  > 

F -<I >=<0PLIS<I>*QPIUS ( I) -2. *Q PLUS (I) +0.75)* A (1*1 1)  - 

C  tCPLUS(I)*GPlUS(IL-a.25)*A  CC2)  . .  .  _ 

10  CONTINUE 
C 

SO  2  G  .  1=  1 .  N-  k. .  _. . 

V  TC  <  D  =  A(I4I1)+F1<:)*R3(I)*  (FOC-ll  +  FTC-l)  )  - 

C  R2(I)*KAC:)*<FO«:)tFTC))  -  2.*f4(:  )*QPLJS(I)  *£PLUSC) 

C  -2.*M  (I)Ml.-QPJLUS  C-l)  )  CFLUSC-l)...  ..  .  _  . 

20  CONTINUE 

VTr  (N-3) =VT0 <N-A> 

o  vT  c  [  o )  =  o ,  c  . .  ....  .  .  _ 

CO  3  0  -  =  1  ,  h-  4 

DVTPtI)r\,TP{;+H-VTGC) 

3C  C  ON T I N  l,J  “  .  ....  .......  . 

o  vT r  ( n- 3 )  =  o.o  • 

c 

CC<  0)  =.  Q.  C  _ _ _ _  ...  _ 

SO  A0  :=1,N-A 
S  =  SIGN  (1  .tCVTCd)  ) 

A1  =  A  MI  M  lS*L¥ia.lIrJJLZR2  (I)/ F.A  (I  )/Zj  ,  APSEDV.ID.il 
C  S*CV’D(I+1)/F2C)/FA(I  )/2,  ) 

FOCI  =  SMMAXllO.tAl) 

A0  CONTINUE . ,  .  .  ... 

C 

cofq:  =  i.n-4 

VC)  =  VTCil)  t  fC(MI’;ltII’C!I)»2,  -  FC  1 1 )  *R2.(  I)  *.R4  1 1  Lf.Z». _ 

50  CONTINUE 

V  (  0  )  =  V  ( 1  ) 

IFC.lQ.U  V  IQ.).  =  0,0...  ...  _ 


.59. 

V  C  N  -  2 )  =  t.  <Ihf'tN-3+ JN) 

40. 

V (N-2)  =  J  (If^  +  N-2  +  JMl 

-1 . 

\l  (N-  1)  =fi  CNM  +  M-1  +  JN) 

42. 

RETURN 

*♦3  . 

£NC 
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